v0.14.0
Loading...
Searching...
No Matches
ComplexConstArea.hpp
Go to the documentation of this file.
1/* This file is part of MoFEM.
2 * MoFEM is free software: you can redistribute it and/or modify it under
3 * the terms of the GNU Lesser General Public License as published by the
4 * Free Software Foundation, either version 3 of the License, or (at your
5 * option) any later version.
6 *
7 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10 * License for more details.
11 *
12 * You should have received a copy of the GNU Lesser General Public
13 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14
15#ifndef __COMPLEX_CONST_AREA_HPP__
16#define __COMPLEX_CONST_AREA_HPP__
17
18#include <math.h>
19#include <complex>
20
21extern "C" {
22
23void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3],
24 double circumcenter[3], double *xi, double *eta,
25 double *zeta);
26void tricircumcenter3d_tp(double a[3], double b[3], double c[3],
27 double circumcenter[3], double *xi, double *eta);
28}
29
31
32/**
33 *
34 * dN/dX = (1/A) * [ Spin[dX/dksi]*dN/deta - Spin[dX/deta]*dN/dksi ]
35 *
36 */
38
40
41 moab::Interface &moab;
42
43 Mat C; //<< constrains matrix
44 Mat Q; //<< projection matrix (usually this is shell matrix)
47
48 ConstrainConstantAarea(MoFEM::Interface &m_field, Mat _C, Mat _Q,
49 string lambda_field_name, int _verbose = 0)
50 : mField(m_field), moab(m_field.get_moab()), C(_C), Q(_Q),
51 lambdaFieldName(lambda_field_name), verbose(_verbose) {
52 // calculate face shape functions direvatives
53 diffNTRI.resize(3, 2);
54 ShapeDiffMBTRI(&*diffNTRI.data().begin());
55 // shape functions Gauss integration weigths
57 // nodal material positions
58 dofs_X.resize(9);
59 // noal values of Lagrange multipliers
60 lambda.resize(3);
61 // dofs indices for Lagrnage multipliers
62 lambdaDofsRowIdx.resize(3);
63 // lambda_dofs_row_ents.resize(3);
64 lambdaDofsColIdx.resize(3);
65 // dofs indices for rows and columns
66 dispDofsRowIdx.resize(9);
67 dispDofsColIdx.resize(9);
68 localDispDofsRowIdx.resize(9);
69 // face node coordinates
70 coords.resize(9);
71 }
72
73 virtual ~ConstrainConstantAarea() = default;
74
75 // elem data
76 MatrixDouble diffNTRI;
77 const double *G_TRI_W;
78 vector<DofIdx> DirichletBC;
82 ublas::vector<double, ublas::bounded_array<double, 9>> coords;
83 ublas::vector<double, ublas::bounded_array<double, 9>> dofs_X;
84 ublas::vector<double, ublas::bounded_array<double, 3>> lambda;
85 // vector<EntityHandle> lambda_dofs_row_ents;
87
88 template <typename DOFS>
89 inline auto loIt(DOFS dofs, const FieldBitNumber bit,
90 const EntityHandle ent) {
91 return dofs->lower_bound(DofEntity::getLoFieldEntityUId(bit, ent));
92 };
93
94 template <typename DOFS>
95 inline auto hiIt(DOFS dofs, const FieldBitNumber bit,
96 const EntityHandle ent) {
97 return dofs->upper_bound(DofEntity::getHiFieldEntityUId(bit, ent));
98 };
99
100 /**
101 * \brief get face data, indices and coors and nodal values
102 *
103 * \param is_that_C_otherwise_dC
104 */
105 PetscErrorCode getData(bool is_that_C_otherwise_dC, bool trans) {
107
108 face = numeredEntFiniteElementPtr->getEnt();
109 fill(lambdaDofsRowIdx.begin(), lambdaDofsRowIdx.end(), -1);
110 fill(lambdaDofsColIdx.begin(), lambdaDofsColIdx.end(), -1);
111 fill(dispDofsRowIdx.begin(), dispDofsRowIdx.end(), -1);
112 fill(dispDofsColIdx.begin(), dispDofsColIdx.end(), -1);
113 fill(localDispDofsRowIdx.begin(), localDispDofsRowIdx.end(), -1);
114 fill(lambda.begin(), lambda.end(), 0);
115 // fill(lambda_dofs_row_ents.begin(),lambda_dofs_row_ents.end(),no_handle);
116 const EntityHandle *conn_face;
117 int num_nodes;
118 CHKERR moab.get_connectivity(face, conn_face, num_nodes, true);
119 if (num_nodes != 3)
120 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
121 "face should have three nodes");
122 CHKERR moab.get_coords(conn_face, num_nodes, &*coords.data().begin());
123
124 const auto field_bit_number = mField.get_field_bit_number(lambdaFieldName);
125 const auto pos_bit_number =
126 mField.get_field_bit_number("MESH_NODE_POSITIONS");
127
128 auto row_dofs = getRowDofsPtr();
129
130 for (int nn = 0; nn < num_nodes; nn++) {
131
132 if (is_that_C_otherwise_dC) {
133
134 auto dit = loIt(row_dofs, field_bit_number, conn_face[nn]);
135 auto hi_dit = hiIt(row_dofs, field_bit_number, conn_face[nn]);
136
137 // it is C
138 // get rows which are Lagrange multipliers
139 if (std::distance(dit, hi_dit) > 0) {
140 if (std::distance(dit, hi_dit) != 1) {
141 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
142 "data inconsistency, number of dof on node for < %s > "
143 "should be 1",
144 lambdaFieldName.c_str());
145 }
146 if (dit->get()->getPetscLocalDofIdx() < 0) {
147 SETERRQ(
148 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
149 "data inconsistency, negative index of local dofs on element");
150 }
151 lambda[nn] = dit->get()->getFieldData();
152 lambdaDofsRowIdx[nn] = dit->get()->getPetscGlobalDofIdx();
153 // lambda_dofs_row_ents[nn] = dit->getEnt();
154 }
155 }
156 if ((!is_that_C_otherwise_dC) || (trans)) {
157
158 auto dit = loIt(row_dofs, field_bit_number, conn_face[nn]);
159 auto hi_dit = hiIt(row_dofs, field_bit_number, conn_face[nn]);
160
161 // it is dC
162 // get rows which are material nodal positions
163 if (std::distance(dit, hi_dit) > 0) {
164 if (std::distance(dit, hi_dit) != 1)
165 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
166 "data inconsistency, number of dof on node for < %s > "
167 "should be 1",
168 lambdaFieldName.c_str());
169 lambda[nn] = dit->get()->getFieldData();
170
171 auto diit = loIt(row_dofs, pos_bit_number, conn_face[nn]);
172 auto hi_diit = hiIt(row_dofs, pos_bit_number, conn_face[nn]);
173
174 if (std::distance(diit, hi_diit) != 3)
175 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
176 "data inconsistency, number of dof on node for "
177 "MESH_NODE_POSITIONS should be 3");
178 for (; diit != hi_diit; diit++) {
179 if (diit->get()->getPetscLocalDofIdx() < 0)
180 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
181 "data inconsistency, negative index of local dofs on "
182 "element");
183 assert(nn * 3 + diit->get()->getDofCoeffIdx() < 9);
184 dispDofsRowIdx[nn * 3 + diit->get()->getDofCoeffIdx()] =
185 diit->get()->getPetscGlobalDofIdx();
186 localDispDofsRowIdx[nn * 3 + diit->get()->getDofCoeffIdx()] =
187 diit->get()->getPetscLocalDofIdx();
188 }
189 }
190 }
191
192
193 auto col_dofs = getColDofsPtr();
194
195 // get columns which are material nodal positions
196 auto dit = loIt(col_dofs, pos_bit_number, conn_face[nn]);
197 auto hi_dit = hiIt(col_dofs, pos_bit_number, conn_face[nn]);
198 if (std::distance(dit, hi_dit) != 3) {
199 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
200 "data inconsistency, number of dof on node for "
201 "MESH_NODE_POSITIONS should be 3, nb dofs = %d",
202 std::distance(dit, hi_dit));
203 }
204 for (; dit != hi_dit; dit++) {
205 if (dit->get()->getPetscLocalDofIdx() < 0) {
206 SETERRQ(
207 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
208 "data inconsistency, negative index of local dofs on element");
209 }
210 dofs_X[nn * 3 + dit->get()->getDofCoeffIdx()] =
211 dit->get()->getFieldData();
212 assert(nn * 3 + dit->get()->getDofCoeffIdx() < 9);
213 dispDofsColIdx[nn * 3 + dit->get()->getDofCoeffIdx()] =
214 dit->get()->getPetscGlobalDofIdx();
215 }
216
217 if (trans) {
218 auto dit = loIt(col_dofs, field_bit_number, conn_face[nn]);
219 auto hi_dit = hiIt(col_dofs, field_bit_number, conn_face[nn]);
220 if (std::distance(dit, hi_dit) > 0) {
221 if (std::distance(dit, hi_dit) != 1) {
222 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
223 "data inconsistency, number of dof on node for < %s > "
224 "should be 1",
225 lambdaFieldName.c_str());
226 }
227 if (dit->get()->getPetscLocalDofIdx() < 0) {
228 SETERRQ(
229 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
230 "data inconsistency, negative index of local dofs on element");
231 }
232 lambdaDofsColIdx[nn] = dit->get()->getPetscGlobalDofIdx();
233 }
234 }
235 }
237 }
238
240
241 /**
242 * \brief calculate direvatives
243 *
244 */
245 PetscErrorCode calcDirevatives(double *diffNTRI, double *dofs_X,
246 double *dofs_iX, double *C, double *iC,
247 double *T, double *iT) {
249 double diffX_xi[3], diffX_eta[3];
250 bzero(diffX_xi, 3 * sizeof(double));
251 bzero(diffX_eta, 3 * sizeof(double));
252 double i_diffX_xi[3], i_diffX_eta[3];
253 bzero(i_diffX_xi, 3 * sizeof(double));
254 bzero(i_diffX_eta, 3 * sizeof(double));
255 Range adj_side_elems;
256
257 BitRefLevel bit = problemPtr->getBitRefLevel();
258 if (!nInTheLoop) {
259 problemTets.clear();
260 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
261 bit, BitRefLevel().set(), MBTET, problemTets);
262 }
263
264 CHKERR mField.getInterface<BitRefManager>()->getAdjacencies(
265 bit, &face, 1, 3, adj_side_elems);
266 adj_side_elems = intersect(problemTets, adj_side_elems);
267 // adj_side_elems = adj_side_elems.subset_by_type(MBTET);
268 if (adj_side_elems.size() == 0) {
269 Range adj_tets_on_surface;
270 BitRefLevel bit_tet_on_surface;
271 bit_tet_on_surface.set(BITREFLEVEL_SIZE - 2);
272 CHKERR mField.getInterface<BitRefManager>()->getAdjacencies(
273 bit_tet_on_surface, &face, 1, 3, adj_tets_on_surface,
274 moab::Interface::INTERSECT, 0);
275
276 adj_side_elems.insert(*adj_tets_on_surface.begin());
277 }
278 if (adj_side_elems.size() != 1) {
279 adj_side_elems.clear();
280 CHKERR mField.getInterface<BitRefManager>()->getAdjacencies(
281 bit, &face, 1, 3, adj_side_elems, moab::Interface::INTERSECT, 5);
282
283 Range::iterator it = adj_side_elems.begin();
284 for (; it != adj_side_elems.end(); it++) {
285 Range nodes;
286 CHKERR mField.get_moab().get_connectivity(&*it, 1, nodes, true);
287 PetscPrintf(PETSC_COMM_WORLD, "%lu %lu %lu %lu\n", nodes[0], nodes[1],
288 nodes[2], nodes[3]);
289 }
290 ParallelComm *pcomm =
291 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
292 if (pcomm->rank() == 0) {
293 EntityHandle out_meshset;
294 CHKERR mField.get_moab().create_meshset(MESHSET_SET, out_meshset);
295 CHKERR mField.get_moab().add_entities(out_meshset, adj_side_elems);
296 CHKERR mField.get_moab().add_entities(out_meshset, &face, 1);
297 CHKERR mField.get_moab().write_file("debug_error.vtk", "VTK", "",
298 &out_meshset, 1);
299 CHKERR mField.get_moab().delete_entities(&out_meshset, 1);
300 }
301 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
302 "expect 1 tet but is %u", adj_side_elems.size());
303 }
304 EntityHandle side_elem = *adj_side_elems.begin();
305 int order[] = {0, 1, 2};
306 if (side_elem != 0) {
307 int side_number, sense, offset;
308 CHKERR moab.side_number(side_elem, face, side_number, sense, offset);
309 if (sense == -1) {
310 order[0] = 1;
311 order[1] = 0;
312 }
313 }
314 // calculate tangent vectors
315 // those vectors are in plane of face
316 for (int nn = 0; nn < 3; nn++) {
317 diffX_xi[0] +=
318 dofs_X[3 * order[nn] + 0] * diffNTRI[2 * nn + 0]; // unit [ m ]
319 diffX_xi[1] += dofs_X[3 * order[nn] + 1] * diffNTRI[2 * nn + 0];
320 diffX_xi[2] += dofs_X[3 * order[nn] + 2] * diffNTRI[2 * nn + 0];
321 diffX_eta[0] += dofs_X[3 * order[nn] + 0] * diffNTRI[2 * nn + 1];
322 diffX_eta[1] += dofs_X[3 * order[nn] + 1] * diffNTRI[2 * nn + 1];
323 diffX_eta[2] += dofs_X[3 * order[nn] + 2] * diffNTRI[2 * nn + 1];
324 if (dofs_iX == NULL)
325 continue;
326 i_diffX_xi[0] += dofs_iX[3 * order[nn] + 0] * diffNTRI[2 * nn + 0];
327 i_diffX_xi[1] += dofs_iX[3 * order[nn] + 1] * diffNTRI[2 * nn + 0];
328 i_diffX_xi[2] += dofs_iX[3 * order[nn] + 2] * diffNTRI[2 * nn + 0];
329 i_diffX_eta[0] += dofs_iX[3 * order[nn] + 0] * diffNTRI[2 * nn + 1];
330 i_diffX_eta[1] += dofs_iX[3 * order[nn] + 1] * diffNTRI[2 * nn + 1];
331 i_diffX_eta[2] += dofs_iX[3 * order[nn] + 2] * diffNTRI[2 * nn + 1];
332 }
333 // spins
334 double SpinX_xi[9], SpinX_eta[9];
335 CHKERR Spin(SpinX_xi, diffX_xi);
336 CHKERR Spin(SpinX_eta, diffX_eta);
337 double iSpinX_xi[9], iSpinX_eta[9];
338 CHKERR Spin(iSpinX_xi, i_diffX_xi);
339 CHKERR Spin(iSpinX_eta, i_diffX_eta);
340 __CLPK_doublecomplex xSpinX_xi[9], xSpinX_eta[9];
341 CHKERR make_complex_matrix(SpinX_xi, iSpinX_xi, xSpinX_xi); // unit [ m ]
342 CHKERR make_complex_matrix(SpinX_eta, iSpinX_eta, xSpinX_eta); // unit [ m ]
343 // calculate complex face normal vector
344 __CLPK_doublecomplex x_one = {1, 0};
345 __CLPK_doublecomplex x_zero = {0, 0};
346 __CLPK_doublecomplex x_normal[3];
347 __CLPK_doublecomplex x_diffX_eta[3];
348 x_diffX_eta[0].r = diffX_eta[0];
349 x_diffX_eta[0].i = i_diffX_eta[0];
350 x_diffX_eta[1].r = diffX_eta[1];
351 x_diffX_eta[1].i = i_diffX_eta[1];
352 x_diffX_eta[2].r = diffX_eta[2];
353 x_diffX_eta[2].i = i_diffX_eta[2];
354 cblas_zgemv(CblasRowMajor, CblasNoTrans, 3, 3, &x_one, xSpinX_xi, 3,
355 x_diffX_eta, 1, &x_zero, x_normal, 1); // unit [ m^2 ]
356 /*//debug
357 Tag th_normal1,th_normal2;
358 double def_NORMAL[] = {0, 0, 0};
359 CHKERR
360 moab.tag_get_handle("NORMAL1",3,MB_TYPE_DOUBLE,th_normal1,MB_TAG_CREAT|MB_TAG_SPARSE,def_NORMAL);
361 CHKERR
362 moab.tag_get_handle("NORMAL2",3,MB_TYPE_DOUBLE,th_normal2,MB_TAG_CREAT|MB_TAG_SPARSE,def_NORMAL);
363 double normal1[3] = { x_normal[0].r, x_normal[1].r, x_normal[2].r };
364 double normal2[3];
365 CHKERR ShapeFaceNormalMBTRI(&diffNTRI[0],&*coords.data().begin(),normal2);
366 if(side_elem!=0) {
367 int side_number,sense,offset;
368 CHKERR moab.side_number(side_elem,face,side_number,sense,offset);
369 if(sense == -1) {
370 cblas_dscal(3,-1,normal2,1);
371 }
372 }
373 double normal1_nrm2 = cblas_dnrm2(3,normal1,1);
374 cblas_dscal(3,1./normal1_nrm2,normal1,1);
375 double normal2_nrm2 = cblas_dnrm2(3,normal2,1);
376 cblas_dscal(3,1./normal2_nrm2,normal2,1);
377 CHKERR moab.tag_set_data(th_normal1,&face,1,normal1);
378 CHKERR moab.tag_set_data(th_normal2,&face,1,normal2); */
379 // calculate complex normal length
380 // double __complex__ x_nrm2 = csqrt(
381 // cpow((x_normal[0].r+I*x_normal[0].i),2+I*0)+
382 // cpow((x_normal[1].r+I*x_normal[1].i),2+I*0)+
383 // cpow((x_normal[2].r+I*x_normal[2].i),2+I*0)
384 // );
385 complex<double> x_nrm2;
386 x_nrm2 = sqrt(pow(complex<double>(x_normal[0].r, x_normal[0].i), 2) +
387 pow(complex<double>(x_normal[1].r, x_normal[1].i), 2) +
388 pow(complex<double>(x_normal[2].r, x_normal[2].i), 2));
389 // calculate dA/dX
390 __CLPK_doublecomplex xNSpinX_xi[3], xNSpinX_eta[3];
391 bzero(xNSpinX_xi, 3 * sizeof(__CLPK_doublecomplex));
392 bzero(xNSpinX_eta, 3 * sizeof(__CLPK_doublecomplex));
393 // __CLPK_doublecomplex x_scalar = { -creal(1./x_nrm2), -cimag(1./x_nrm2) };
394 // // unit [ 1/m^2 ]
395 __CLPK_doublecomplex x_scalar = {-real(1. / x_nrm2),
396 -imag(1. / x_nrm2)}; // unit [ 1/m^2 ]
397 cblas_zgemv(CblasRowMajor, CblasNoTrans, 3, 3, &x_scalar, xSpinX_xi, 3,
398 x_normal, 1, &x_zero, xNSpinX_xi,
399 1); // unit [ 1/m^2 * m = 1/m ]
400 cblas_zgemv(CblasRowMajor, CblasNoTrans, 3, 3, &x_scalar, xSpinX_eta, 3,
401 x_normal, 1, &x_zero, xNSpinX_eta,
402 1); // unit [ 1/m^2 * m = 1/m ]
403 if (C != NULL)
404 bzero(C, 9 * sizeof(double));
405 if (iC != NULL)
406 bzero(iC, 9 * sizeof(double));
407 if (T != NULL)
408 bzero(T, 9 * sizeof(double));
409 if (iT != NULL)
410 bzero(iT, 9 * sizeof(double));
411 for (int nn = 0; nn < 3; nn++) {
412 double A[3], iA[3];
413 A[0] = xNSpinX_xi[0].r * diffNTRI[2 * order[nn] + 1] -
414 xNSpinX_eta[0].r * diffNTRI[2 * order[nn] + 0]; // unit [ 1/m ]
415 A[1] = xNSpinX_xi[1].r * diffNTRI[2 * order[nn] + 1] -
416 xNSpinX_eta[1].r * diffNTRI[2 * order[nn] + 0];
417 A[2] = xNSpinX_xi[2].r * diffNTRI[2 * order[nn] + 1] -
418 xNSpinX_eta[2].r * diffNTRI[2 * order[nn] + 0];
419 iA[0] = xNSpinX_xi[0].i * diffNTRI[2 * order[nn] + 1] -
420 xNSpinX_eta[0].i * diffNTRI[2 * order[nn] + 0]; // unit [ 1/m ]
421 iA[1] = xNSpinX_xi[1].i * diffNTRI[2 * order[nn] + 1] -
422 xNSpinX_eta[1].i * diffNTRI[2 * order[nn] + 0];
423 iA[2] = xNSpinX_xi[2].i * diffNTRI[2 * order[nn] + 1] -
424 xNSpinX_eta[2].i * diffNTRI[2 * order[nn] + 0];
425 if (C != NULL) {
426 C[3 * nn + 0] = A[0];
427 C[3 * nn + 1] = A[1];
428 C[3 * nn + 2] = A[2];
429 }
430 if (iC != NULL) {
431 iC[3 * nn + 0] = iA[0];
432 iC[3 * nn + 1] = iA[1];
433 iC[3 * nn + 2] = iA[2];
434 }
435 if ((T != NULL) || (iT != NULL)) {
436 double SpinA[9];
437 CHKERR Spin(SpinA, A); // unit [1/m]
438 double iSpinA[9];
439 CHKERR Spin(iSpinA, iA); // unit [1/m]
440 __CLPK_doublecomplex xSpinA[9];
441 // make spin matrix to calculate cross product
442 CHKERR make_complex_matrix(SpinA, iSpinA, xSpinA);
444 cblas_zgemv(CblasRowMajor, CblasNoTrans, 3, 3, &x_scalar, xSpinA, 3,
445 x_normal, 1, &x_zero, xT, 1); // unit [1/m]
446 if (T != NULL) {
447 T[3 * nn + 0] = xT[0].r;
448 T[3 * nn + 1] = xT[1].r;
449 T[3 * nn + 2] = xT[2].r;
450 }
451 if (iT != NULL) {
452 iT[3 * nn + 0] = xT[0].i;
453 iT[3 * nn + 1] = xT[1].i;
454 iT[3 * nn + 2] = xT[2].i;
455 }
456 }
457 }
458 if (C != NULL)
459 cblas_dscal(9, 0.25, C, 1);
460 if (iC != NULL)
461 cblas_dscal(9, 0.25, iC, 1);
462 if (T != NULL)
463 cblas_dscal(9, 0.25, T, 1);
464 if (iT != NULL)
465 cblas_dscal(9, 0.25, iT, 1);
467 }
468
469 map<DofIdx, Vec> mapV;
470 PetscErrorCode preProcess() {
472
473 if (Q != PETSC_NULL) {
474 const auto field_bit_number =
477 problemPtr, field_bit_number, dofs_it)) {
478 CHKERR MatCreateVecs(C, &mapV[dofs_it->get()->getPetscGlobalDofIdx()],
479 PETSC_NULL);
480 CHKERR VecZeroEntries(mapV[dofs_it->get()->getPetscGlobalDofIdx()]);
481 }
482 }
483
485 }
486
487 PetscErrorCode operator()() {
489
490 CHKERR getData(true, false);
491
492 ublas::vector<double, ublas::bounded_array<double, 9>> ELEM_CONSTRAIN(9);
493 CHKERR calcDirevatives(&*diffNTRI.data().begin(), &*dofs_X.data().begin(),
494 NULL, &*ELEM_CONSTRAIN.data().begin(), NULL, NULL,
495 NULL);
496 for (int nn = 0; nn < 3; nn++) {
497 if (lambdaDofsRowIdx[nn] == -1)
498 continue;
499 for (int NN = 0; NN < 3; NN++) {
500 if (NN != nn)
501 continue;
502 if (Q == PETSC_NULL) {
503 CHKERR MatSetValues(C, 1, &(lambdaDofsRowIdx.data()[nn]), 3,
504 &(dispDofsColIdx.data()[3 * NN]),
505 &ELEM_CONSTRAIN.data()[3 * NN], ADD_VALUES);
506 }
507 if (Q != PETSC_NULL) {
508 if (mapV.find(lambdaDofsRowIdx[nn]) == mapV.end())
509 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
510 CHKERR VecSetValues(mapV[lambdaDofsRowIdx[nn]], 3,
511 &(dispDofsColIdx.data()[3 * NN]),
512 &ELEM_CONSTRAIN.data()[3 * NN], ADD_VALUES);
513 }
514 }
515 }
516
518 }
519
520 PetscErrorCode postProcess() {
522 if (Q != PETSC_NULL) {
523 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
524 Vec Qv;
525 CHKERR mField.getInterface<VecManager>()->vecCreateGhost(
526 problemPtr->getName(), COL, &Qv);
527 const auto field_bit_number =
530 field_bit_number, dofs)) {
531 if (mapV.find(dofs->get()->getPetscGlobalDofIdx()) == mapV.end()) {
532 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
533 "data inconsistency");
534 }
535 CHKERR VecAssemblyBegin(mapV[dofs->get()->getPetscGlobalDofIdx()]);
536 CHKERR VecAssemblyEnd(mapV[dofs->get()->getPetscGlobalDofIdx()]);
537 CHKERR MatMult(Q, mapV[dofs->get()->getPetscGlobalDofIdx()], Qv);
538 CHKERR VecGhostUpdateBegin(Qv, INSERT_VALUES, SCATTER_FORWARD);
539 CHKERR VecGhostUpdateEnd(Qv, INSERT_VALUES, SCATTER_FORWARD);
540 double *array;
541 CHKERR VecGetArray(Qv, &array);
542 vector<DofIdx> glob_idx;
543 vector<double> vals;
545 dofs->get()->getEnt(), diit)) {
546 if (diit->get()->getPart() != pcomm->rank())
547 continue;
548 if (diit->get()->getName() != "MESH_NODE_POSITIONS")
549 continue;
550 glob_idx.push_back(diit->get()->getPetscGlobalDofIdx());
551 if (diit->get()->getPetscLocalDofIdx() < 0) {
552 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
553 "data inconsistency");
554 }
555 vals.push_back(array[diit->get()->getPetscLocalDofIdx()]);
556 }
557 int row = dofs->get()->getPetscGlobalDofIdx();
558 CHKERR MatSetValues(C, 1, &row, glob_idx.size(),
559 glob_idx.empty() ? PETSC_NULL : &glob_idx[0],
560 vals.empty() ? PETSC_NULL : &vals[0], ADD_VALUES);
561 CHKERR VecRestoreArray(Qv, &array);
562 CHKERR VecDestroy(&mapV[dofs->get()->getPetscGlobalDofIdx()]);
563 }
564 CHKERR VecDestroy(&Qv);
565 mapV.clear();
566 }
568 }
569};
570
571/**
572 * Constrains is vector_tangent_to_front*nodal_element_quality_at_crack_front =
573 * 0, multiplying that constrain by Legrange multipliers, we get
574 *
575 * 1) vector_tangent_to_front*nodal_element_quality_at_crack_front=0 for each
576 * lambda, 2) lambda*vector_tangent_to_front*direvatives_of_quality = 0
577 *
578 * Those two equations need to be linearsed and assembled into system of
579 * equations
580 *
581 *
582 */
584
587 const double eps;
589
591 string lambda_field_name,
592 int verbose = 0)
593 : ConstrainConstantAarea(m_field, PETSC_NULL, PETSC_NULL,
594 lambda_field_name, verbose),
595 frontF(PETSC_NULL), tangentFrontF(PETSC_NULL), eps(1e-10),
596 ownVectors(false) {}
597
599 if (ownVectors) {
600 CHKERR VecDestroy(&frontF);
601 CHKERRABORT(PETSC_COMM_WORLD, ierr);
602 CHKERR VecDestroy(&tangentFrontF);
603 CHKERRABORT(PETSC_COMM_WORLD, ierr);
604 }
605 }
606
608 PetscErrorCode preProcess() {
610
611 switch (ts_ctx) {
612 case CTX_TSSETIFUNCTION: {
613 snes_ctx = CTX_SNESSETFUNCTION;
614 snes_f = ts_F;
615 break;
616 }
617 case CTX_TSSETIJACOBIAN: {
618 snes_ctx = CTX_SNESSETJACOBIAN;
619 snes_B = ts_B;
620 break;
621 }
622 default:
623 break;
624 }
625
627
628 switch (snes_ctx) {
629 case CTX_SNESSETFUNCTION: {
630 CHKERR VecAssemblyBegin(snes_f);
631 CHKERR VecAssemblyEnd(snes_f);
632 CHKERR VecZeroEntries(tangentFrontF);
633 CHKERR VecGhostUpdateBegin(tangentFrontF, INSERT_VALUES, SCATTER_FORWARD);
634 CHKERR VecGhostUpdateEnd(tangentFrontF, INSERT_VALUES, SCATTER_FORWARD);
635 CHKERR VecAssemblyBegin(tangentFrontF);
636 CHKERR VecAssemblyEnd(tangentFrontF);
637 } break;
638 case CTX_SNESSETJACOBIAN: {
639 CHKERR MatAssemblyBegin(snes_B, MAT_FLUSH_ASSEMBLY);
640 CHKERR MatAssemblyEnd(snes_B, MAT_FLUSH_ASSEMBLY);
641 } break;
642 default: {
643 } break;
644 }
646 }
647
648 PetscErrorCode operator()() {
650 // get elemet hanlde
651 EntityHandle face = numeredEntFiniteElementPtr->getEnt();
652 // get adjacent tets
653 Range tet;
654 CHKERR mField.get_moab().get_adjacencies(&face, 1, 3, false, tet);
655 // get data and indices on dofs
656 CHKERR getData(true, true);
657
658 // Calculate tangent vector constrains
659 ublas::vector<double, ublas::bounded_array<double, 9>> tangent_constrains(
660 9);
661 CHKERR calcDirevatives(&*diffNTRI.data().begin(), &*dofs_X.data().begin(),
662 NULL, NULL, NULL,
663 &*tangent_constrains.data().begin(), NULL);
664 // take in account face orientation in respect crack surface
665 Tag th_interface_side;
666 CHKERR moab.tag_get_handle("INTERFACE_SIDE", th_interface_side);
667 int side;
668 CHKERR moab.tag_get_data(th_interface_side, &face, 1, &side);
669 if (side == 1) {
670 tangent_constrains *= -1;
671 }
672 // get smoothing element forces
673 ublas::vector<double, ublas::bounded_array<double, 9>>
674 f_front_mesh_smoothing(9);
675 ublas::noalias(f_front_mesh_smoothing) = ublas::zero_vector<double>(9);
676 double *f_front_mesh_array;
677 CHKERR VecGetArray(frontF, &f_front_mesh_array);
678 for (int nn = 0; nn < 3; nn++) {
679 for (int dd = 0; dd < 3; dd++) {
680 if (localDispDofsRowIdx[3 * nn + dd] == -1)
681 continue;
682 f_front_mesh_smoothing[3 * nn + dd] =
683 f_front_mesh_array[localDispDofsRowIdx[3 * nn + dd]];
684 }
685 }
686 CHKERR VecRestoreArray(frontF, &f_front_mesh_array);
687 // calculate tangent matrix
688 if (snes_ctx == CTX_SNESSETJACOBIAN) {
689 // get radius
690 double center[3];
691 tricircumcenter3d_tp(&coords.data()[0], &coords.data()[3],
692 &coords.data()[6], center, NULL, NULL);
693 cblas_daxpy(3, -1, &coords.data()[0], 1, center, 1);
694 double r = cblas_dnrm2(3, center, 1);
695 for (int nn = 0; nn < 3; nn++) {
696 for (int dd = 0; dd < 3; dd++) {
697 // ---> calculate tangent starts here
698 ublas::vector<double, ublas::bounded_array<double, 9>> idofs_X(9, 0);
699 idofs_X[nn * 3 + dd] = r * eps;
700 ublas::vector<double, ublas::bounded_array<double, 9>>
701 direvative_of_constrains(9);
702 // use complex derivative
703 CHKERR calcDirevatives(&*diffNTRI.data().begin(),
704 &*dofs_X.data().begin(),
705 &*idofs_X.data().begin(), NULL, NULL, NULL,
706 &*direvative_of_constrains.data().begin());
707 if (side == 1) {
708 direvative_of_constrains /= -r * eps;
709 } else {
710 direvative_of_constrains /= +r * eps;
711 }
712 // dg -> C*q_quality
713 double g[3] = {0, 0, 0};
714 for (int nnn = 0; nnn < 3; nnn++) {
715 if (lambdaDofsRowIdx[nnn] == -1)
716 continue;
717 for (int ddd = 0; ddd < 3; ddd++) {
718 g[nnn] += direvative_of_constrains[nnn * 3 + ddd] *
719 f_front_mesh_smoothing[3 * nnn + ddd];
720 }
721 }
722 for (int nnn = 0; nnn < 3; nnn++) {
723 if (lambdaDofsRowIdx[nnn] == -1)
724 continue;
725 CHKERR MatSetValues(snes_B, 1, &lambdaDofsRowIdx[nnn], 1,
726 &dispDofsColIdx[3 * nn + dd], &g[nnn],
727 ADD_VALUES);
728 }
729 // dCT_lambda
730 for (int nnn = 0; nnn < 3; nnn++) {
731 for (int ddd = 0; ddd < 3; ddd++) {
732 direvative_of_constrains[nnn * 3 + ddd] *= lambda[nnn];
733 }
734 }
735 for (int nnn = 0; nnn < 3; nnn++) {
736 if (lambdaDofsRowIdx[nnn] == -1)
737 continue;
738 CHKERR MatSetValues(snes_B, 3, &dispDofsRowIdx[3 * nnn], 1,
739 &dispDofsColIdx[3 * nn + dd],
740 &direvative_of_constrains[3 * nnn], ADD_VALUES);
741 }
742 // ---> calculate tangent end here
743 }
744 }
745 }
746 switch (snes_ctx) {
747 case CTX_SNESSETFUNCTION: {
748 ublas::vector<double, ublas::bounded_array<double, 3>> g(3);
749 for (int nn = 0; nn < 3; nn++) {
750 g[nn] = cblas_ddot(3, &tangent_constrains[3 * nn], 1,
751 &f_front_mesh_smoothing[3 * nn], 1);
752 }
753 // cerr << "g : " << g << endl;
754 CHKERR VecSetValues(snes_f, 3, &*lambdaDofsRowIdx.data().begin(),
755 &*g.data().begin(), ADD_VALUES);
756 CHKERR VecSetValues(tangentFrontF, 9, &dispDofsRowIdx[0],
757 &*tangent_constrains.data().begin(), ADD_VALUES);
758 ublas::vector<double, ublas::bounded_array<double, 9>> f(9);
759 for (int nn = 0; nn < 3; nn++) {
760 for (int dd = 0; dd < 3; dd++) {
761 f[nn * 3 + dd] = lambda[nn] * tangent_constrains[3 * nn + dd];
762 }
763 }
764 // cerr << "f : " << f << endl;
765 CHKERR VecSetValues(snes_f, 9, &dispDofsRowIdx[0], &*f.data().begin(),
766 ADD_VALUES);
767 /*//TAG - only for one proc analysis
768 for(int nn = 0;nn<3;nn++) {
769 EntityHandle ent = lambda_dofs_row_ents[nn];
770 if(ent == no_handle) continue;
771 double *t;
772 CHKERR mField.get_moab().tag_get_by_ptr(thFrontTangent,&ent,1,(const void
773 **)&t); cblas_daxpy(3,+1,&tangent_constrains[3*nn],1,t,1);
774 }*/
775 } break;
776 case CTX_SNESSETJACOBIAN: {
777 for (int nn = 0; nn < 3; nn++) {
778 int lambda_dof_idx = lambdaDofsColIdx[nn];
779 CHKERR MatSetValues(snes_B, 3, &dispDofsRowIdx[3 * nn], 1,
780 &lambda_dof_idx, &tangent_constrains[3 * nn],
781 ADD_VALUES);
782 }
783 } break;
784 default:
785 break;
786 }
788 }
789
790 PetscErrorCode postProcess() {
792 switch (snes_ctx) {
793 case CTX_SNESSETFUNCTION: {
794 CHKERR VecAssemblyBegin(tangentFrontF);
795 CHKERR VecAssemblyEnd(tangentFrontF);
796 CHKERR VecGhostUpdateBegin(tangentFrontF, ADD_VALUES, SCATTER_REVERSE);
797 CHKERR VecGhostUpdateEnd(tangentFrontF, ADD_VALUES, SCATTER_REVERSE);
798 CHKERR VecGhostUpdateBegin(tangentFrontF, INSERT_VALUES, SCATTER_FORWARD);
799 CHKERR VecGhostUpdateEnd(tangentFrontF, INSERT_VALUES, SCATTER_FORWARD);
800 } break;
801 case CTX_SNESSETJACOBIAN: {
802 CHKERR MatAssemblyBegin(snes_B, MAT_FLUSH_ASSEMBLY);
803 CHKERR MatAssemblyEnd(snes_B, MAT_FLUSH_ASSEMBLY);
804 } break;
805 default: {
806 } break;
807 }
809 }
810};
811
812} // namespace ObosleteUsersModules
813
814#endif // __COMPLEX_CONST_AREA_HPP__
void tricircumcenter3d_tp(double a[3], double b[3], double c[3], double circumcenter[3], double *xi, double *eta)
void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3], double circumcenter[3], double *xi, double *eta, double *zeta)
static PetscErrorCode ierr
@ COL
Definition: definitions.h:123
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:219
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:194
static const double G_TRI_W1[]
Definition: fem_tools.h:314
PetscErrorCode Spin(double *spinOmega, double *vecOmega)
calculate spin matrix from vector
Definition: fem_tools.c:546
PetscErrorCode make_complex_matrix(double *reA, double *imA, __CLPK_doublecomplex *xA)
Compose complex matrix (3x3) from two real matrices.
Definition: fem_tools.c:557
double eta
#define _IT_NUMEREDDOF_ROW_BY_BITNUMBER_FOR_LOOP_(PROBLEMPTR, FIELD_BIT_NUMBER, IT)
#define _IT_NUMEREDDOF_COL_BY_ENT_FOR_LOOP_(PROBLEMPTR, ENT, IT)
use with loops to iterate col DOFs
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
const double T
constexpr AssemblyType A
double zeta
Viscous hardening.
Definition: plastic.cpp:177
constexpr double g
__CLPK_doublereal i
Definition: lapack_wrap.h:35
__CLPK_doublereal r
Definition: lapack_wrap.h:35
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
PetscErrorCode getData(bool is_that_C_otherwise_dC, bool trans)
get face data, indices and coors and nodal values
ublas::vector< double, ublas::bounded_array< double, 9 > > dofs_X
PetscErrorCode calcDirevatives(double *diffNTRI, double *dofs_X, double *dofs_iX, double *C, double *iC, double *T, double *iT)
calculate direvatives
ublas::vector< double, ublas::bounded_array< double, 3 > > lambda
ConstrainConstantAarea(MoFEM::Interface &m_field, Mat _C, Mat _Q, string lambda_field_name, int _verbose=0)
auto hiIt(DOFS dofs, const FieldBitNumber bit, const EntityHandle ent)
auto loIt(DOFS dofs, const FieldBitNumber bit, const EntityHandle ent)
ublas::vector< double, ublas::bounded_array< double, 9 > > coords
TangentWithMeshSmoothingFrontConstrain(MoFEM::Interface &m_field, string lambda_field_name, int verbose=0)
void tricircumcenter3d_tp(a, b, c, circumcenter, double *xi, double *eta)