v0.15.0
Loading...
Searching...
No Matches
SurfacePressureComplexForLazy.cpp
Go to the documentation of this file.
1/* \file SurfacePressureComplexForLazy.cpp
2 */
3
4
5
6#include <MoFEM.hpp>
7
8using namespace MoFEM;
9#include <MethodForForceScaling.hpp>
10#include <SurfacePressure.hpp>
11#include <SurfacePressureComplexForLazy.hpp>
12
13extern "C" {
14
15MoFEMErrorCode Traction_hierarchical(int order, int *order_edge, double *N,
16 double *N_face, double *N_edge[],
17 double *t, double *t_edge[],
18 double *t_face, double *traction, int gg);
19MoFEMErrorCode Fext_h_hierarchical(
20 int order, int *order_edge, double *N, double *N_face, double *N_edge[],
21 double *diffN, double *diffN_face, double *diffN_edge[], double *t,
22 double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[],
23 double *dofs_x_face, double *idofs_x, double *idofs_x_edge[],
24 double *idofs_x_face, double *Fext, double *Fext_egde[], double *Fext_face,
25 double *iFext, double *iFext_egde[], double *iFext_face, int g_dim,
26 const double *g_w);
27MoFEMErrorCode KExt_hh_hierarchical(
28 double eps, int order, int *order_edge, double *N, double *N_face,
29 double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[],
30 double *t, double *t_edge[], double *t_face, double *dofs_x,
31 double *dofs_x_edge[], double *dofs_x_face, double *KExt_hh,
32 double *KExt_egdeh[3], double *KExt_faceh, int g_dim, const double *g_w);
33MoFEMErrorCode KExt_hh_hierarchical_edge(
34 double eps, int order, int *order_edge, double *N, double *N_face,
35 double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[],
36 double *t, double *t_edge[], double *t_face, double *dofs_x,
37 double *dofs_x_edge[], double *dofs_x_face, double *Khext_edge[3],
38 double *KExt_edgeegde[3][3], double *KExt_faceedge[3], int g_dim,
39 const double *g_w);
40MoFEMErrorCode
41KExt_hh_hierarchical_face(double eps, int order, int *order_edge, double *N,
42 double *N_face, double *N_edge[], double *diffN,
43 double *diffN_face, double *diffN_edge[], double *t,
44 double *t_edge[], double *t_face, double *dofs_x,
45 double *dofs_x_edge[], double *dofs_x_face,
46 double *KExt_hface, double *KExt_egdeface[3],
47 double *KExt_faceface, int g_dim, const double *g_w);
48
49void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3],
50 double circumcenter[3], double *xi, double *eta,
51 double *zeta);
52void tricircumcenter3d_tp(double a[3], double b[3], double c[3],
53 double circumcenter[3], double *xi, double *eta);
54}
55
60
65
66MoFEMErrorCode NeumannForcesSurfaceComplexForLazy::AuxMethodSpatial::doWork(
67 int side, EntityType type, EntitiesFieldData::EntData &data) {
68 MoFEMFunctionBegin;
69
70 switch (type) {
71 case MBVERTEX: {
72 if (data.getFieldData().size() != 9) {
73 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
74 "it should be 9 dofs on vertices but is %ld of field < %s >",
75 data.getFieldData().size(), rowFieldName.c_str());
76 }
77 myPtr->N = &*data.getN().data().begin();
78 myPtr->diffN = &*data.getDiffN().data().begin();
79 myPtr->dOfs_x.resize(data.getFieldData().size());
80 ublas::noalias(myPtr->dOfs_x) = data.getFieldData();
81 myPtr->dofs_x = &*myPtr->dOfs_x.data().begin();
82 myPtr->dOfs_x_indices.resize(data.getIndices().size());
83 ublas::noalias(myPtr->dOfs_x_indices) = data.getIndices();
84 myPtr->dofs_x_indices = &*myPtr->dOfs_x_indices.data().begin();
85 } break;
86 case MBEDGE: {
87 myPtr->order_edge[side] = data.getOrder();
88 myPtr->N_edge[side] = &*data.getN().data().begin();
89 myPtr->diffN_edge[side] = &*data.getDiffN().data().begin();
90 myPtr->dOfs_x_edge.resize(3);
91 myPtr->dOfs_x_edge[side].resize(data.getFieldData().size());
92 myPtr->dofs_x_edge[side] = &*myPtr->dOfs_x_edge[side].data().begin();
93 myPtr->dOfs_x_edge_indices.resize(3);
94 myPtr->dOfs_x_edge_indices[side].resize(data.getIndices().size());
95 ublas::noalias(myPtr->dOfs_x_edge_indices[side]) = data.getIndices();
96 myPtr->dofs_x_edge_indices[side] =
97 &*myPtr->dOfs_x_edge_indices[side].data().begin();
98 } break;
99 case MBTRI: {
100 myPtr->order_face = data.getOrder();
101 myPtr->N_face = &*data.getN().data().begin();
102 myPtr->diffN_face = &*data.getDiffN().data().begin();
103 myPtr->dOfs_x_face.resize(data.getFieldData().size());
104 ublas::noalias(myPtr->dOfs_x_face) = data.getFieldData();
105 myPtr->dofs_x_face = &*myPtr->dOfs_x_face.data().begin();
106 myPtr->dOfs_x_face_indices.resize(data.getIndices().size());
107 ublas::noalias(myPtr->dOfs_x_face_indices) = data.getIndices();
108 myPtr->dofs_x_face_indices = &*myPtr->dOfs_x_face_indices.data().begin();
109 } break;
110 default:
111 SETERRQ(PETSC_COMM_SELF, 1, "unknown entity type");
112 }
113
114 MoFEMFunctionReturn(0);
115}
116
117MoFEMErrorCode NeumannForcesSurfaceComplexForLazy::AuxMethodMaterial::doWork(
118 int side, EntityType type, EntitiesFieldData::EntData &data) {
119 MoFEMFunctionBegin;
120
121 switch (type) {
122 case MBVERTEX: {
123 if (data.getFieldData().size() != 9) {
124 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
125 "it should be 9 dofs on vertices but is %ld",
126 data.getFieldData().size());
127 }
128 if (data.getN().size2() != 3) {
129 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
130 "it should 3 shape functions for 3 nodes");
131 }
132 myPtr->N = &*data.getN().data().begin();
133 myPtr->diffN = &*data.getDiffN().data().begin();
134 myPtr->dOfs_X_indices.resize(data.getIndices().size());
135 ublas::noalias(myPtr->dOfs_X_indices) = data.getIndices();
136 myPtr->dofs_X_indices = &*myPtr->dOfs_X_indices.data().begin();
137 myPtr->dOfs_X.resize(data.getFieldData().size());
138 ublas::noalias(myPtr->dOfs_X) = data.getFieldData();
139 myPtr->dofs_X = &*myPtr->dOfs_X.data().begin();
140 } break;
141 case MBEDGE: {
142 myPtr->order_edge_material[side] = data.getOrder();
143 myPtr->dOfs_X_edge.resize(3);
144 myPtr->dOfs_X_edge[side].resize(data.getFieldData().size());
145 ublas::noalias(myPtr->dOfs_X_edge[side]) = data.getFieldData();
146 myPtr->dofs_X_edge[side] = &*myPtr->dOfs_X_edge[side].data().begin();
147 } break;
148 case MBTRI: {
149 myPtr->order_face_material = data.getOrder();
150 myPtr->dOfs_X_face.resize(data.getFieldData().size());
151 ublas::noalias(myPtr->dOfs_X_face) = data.getFieldData();
152 myPtr->dofs_X_face = &*myPtr->dOfs_X_face.data().begin();
153 } break;
154 default:
155 SETERRQ(PETSC_COMM_SELF, 1, "unknown entity type");
156 }
157
158 MoFEMFunctionReturn(0);
159}
160
161NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::MyTriangleSpatialFE(
162 MoFEM::Interface &_mField, Mat _Aij, Vec &_F, double *scale_lhs,
163 double *scale_rhs, std::string spatial_field_name,
164 std::string mat_field_name)
165 : FaceElementForcesAndSourcesCore(_mField), sCaleLhs(scale_lhs),
166 sCaleRhs(scale_rhs), typeOfForces(CONSERVATIVE), eps(1e-8), uSeF(false) {
167
168 meshPositionsFieldName = "NoNE";
169 methodsOp.clear();
170
171 Aij = _Aij;
172 F = _F;
173
174 snes_B = _Aij;
175 snes_f = _F;
176
177 if (mField.check_field(mat_field_name)) {
178 getOpPtrVector().push_back(new AuxMethodMaterial(
179 mat_field_name, this, ForcesAndSourcesCore::UserDataOperator::OPROW));
180 }
181 getOpPtrVector().push_back(new AuxMethodSpatial(
182 spatial_field_name, this, ForcesAndSourcesCore::UserDataOperator::OPROW));
183}
184
185MoFEMErrorCode NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::rHs() {
186 MoFEMFunctionBegin;
187
188 auto &dataH1 = *dataOnElement[H1];
189
190 fExtNode.resize(9);
191 fExtFace.resize(dataH1.dataOnEntities[MBTRI][0].getFieldData().size());
192 fExtEdge.resize(3);
193 for (int ee = 0; ee < 3; ee++) {
194 int nb_edge_dofs = dOfs_x_edge_indices[ee].size();
195 if (nb_edge_dofs > 0) {
196 fExtEdge[ee].resize(nb_edge_dofs);
197 Fext_edge[ee] = &*fExtEdge[ee].data().begin();
198 } else {
199 Fext_edge[ee] = NULL;
200 }
201 }
202
203 switch (typeOfForces) {
204 case CONSERVATIVE:
205 CHKERR Fext_h_hierarchical(
206 order_face, order_edge, // 2
207 N, N_face, N_edge, diffN, diffN_face, diffN_edge, // 8
208 t_loc, NULL, NULL, // 11
209 dofs_x, dofs_x_edge, dofs_x_face, // 14
210 NULL, NULL, NULL, // 17
211 &*fExtNode.data().begin(), Fext_edge, &*fExtFace.data().begin(), // 20
212 NULL, NULL, NULL, // 23
213 gaussPts.size2(), &gaussPts(2, 0));
214 break;
215 case NONCONSERVATIVE:
216 for (int ee = 0; ee < 3; ee++) {
217 dOfs_X_edge.resize(3);
218 unsigned int s = dOfs_X_edge[ee].size();
219 dOfs_X_edge[ee].resize(dOfs_x_edge[ee].size(), true);
220 for (; s < dOfs_X_edge[ee].size(); s++) {
221 dOfs_X_edge[ee][s] = 0;
222 }
223 dofs_X_edge[ee] = &*dOfs_X_edge[ee].data().begin();
224 }
225 unsigned int s = dOfs_X_face.size();
226 dOfs_X_face.resize(dOfs_x_face.size(), true);
227 for (; s < dOfs_X_face.size(); s++) {
228 dOfs_X_face[s] = 0;
229 }
230 dofs_X_face = &*dOfs_X_face.data().begin();
231
232 CHKERR Fext_h_hierarchical(
233 order_face, order_edge, // 2
234 N, N_face, N_edge, diffN, diffN_face, diffN_edge, // 8
235 t_loc, NULL, NULL, // 11
236 dofs_X, dofs_X_edge, dofs_X_face, // 14
237 NULL, NULL, NULL, // 17
238 &*fExtNode.data().begin(), Fext_edge, &*fExtFace.data().begin(), // 20
239 NULL, NULL, NULL, // 23
240 gaussPts.size2(), &gaussPts(2, 0));
241 break;
242 }
243
244 Vec f = snes_f;
245 if (uSeF)
246 f = F;
247
248 CHKERR VecSetValues(f, 9, dofs_x_indices, &*fExtNode.data().begin(),
249 ADD_VALUES);
250 if (dOfs_x_face_indices.size() > 0) {
251 CHKERR VecSetValues(f, dOfs_x_face_indices.size(), dofs_x_face_indices,
252 &*fExtFace.data().begin(), ADD_VALUES);
253 }
254 for (int ee = 0; ee < 3; ee++) {
255 if (dOfs_x_edge_indices[ee].size() > 0) {
256 CHKERR VecSetValues(f, dOfs_x_edge_indices[ee].size(),
257 dofs_x_edge_indices[ee], Fext_edge[ee], ADD_VALUES);
258 }
259 }
260
261 MoFEMFunctionReturn(0);
262}
263
264MoFEMErrorCode NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::lHs() {
265 MoFEMFunctionBegin;
266
267 if (typeOfForces == NONCONSERVATIVE) {
268 MoFEMFunctionReturnHot(0);
269 }
270
271 auto &dataH1 = *dataOnElement[H1];
272
273 double center[3];
274 tricircumcenter3d_tp(&coords.data()[0], &coords.data()[3], &coords.data()[6],
275 center, NULL, NULL);
276 cblas_daxpy(3, -1, &coords.data()[0], 1, center, 1);
277 double r = cblas_dnrm2(3, center, 1);
278
279 kExtNodeNode.resize(9, 9);
280 kExtEdgeNode.resize(3);
281 for (int ee = 0; ee < 3; ee++) {
282 kExtEdgeNode[ee].resize(dOfs_x_edge_indices[ee].size(), 9);
283 Kext_edge_node[ee] = &*kExtEdgeNode[ee].data().begin();
284 }
285 kExtFaceNode.resize(dOfs_x_face_indices.size(), 9);
286 CHKERR KExt_hh_hierarchical(
287 r * eps, order_face, order_edge, N, N_face, N_edge, diffN, diffN_face,
288 diffN_edge, t_loc, NULL, NULL, dofs_x, dofs_x_edge, dofs_x_face,
289 &*kExtNodeNode.data().begin(), Kext_edge_node,
290 &*kExtFaceNode.data().begin(), gaussPts.size2(), &gaussPts(2, 0));
291 // cerr << kExtNodeNode << endl;
292 CHKERR MatSetValues(snes_B, 9, dofs_x_indices, 9, dofs_x_indices,
293 &*kExtNodeNode.data().begin(), ADD_VALUES);
294 CHKERR MatSetValues(snes_B, kExtFaceNode.size1(), dofs_x_face_indices, 9,
295 dofs_x_indices, &*kExtFaceNode.data().begin(),
296 ADD_VALUES);
297 // cerr << kExtFaceNode << endl;
298 for (int ee = 0; ee < 3; ee++) {
299 // cerr << kExtEdgeNode[ee] << endl;
300 CHKERR MatSetValues(snes_B, kExtEdgeNode[ee].size1(),
301 dofs_x_edge_indices[ee], 9, dofs_x_indices,
302 Kext_edge_node[ee], ADD_VALUES);
303 }
304
305 kExtNodeFace.resize(9, dOfs_x_face_indices.size());
306 kExtEdgeFace.resize(3);
307 for (int ee = 0; ee < 3; ee++) {
308 kExtEdgeFace[ee].resize(
309 dOfs_x_edge_indices[ee].size(),
310 dataH1.dataOnEntities[MBTRI][0].getIndices().size());
311 Kext_edge_face[ee] = &*kExtEdgeFace[ee].data().begin();
312 }
313 kExtFaceFace.resize(dOfs_x_face_indices.size(), dOfs_x_face_indices.size());
314 CHKERR KExt_hh_hierarchical_face(
315 r * eps, order_face, order_edge, N, N_face, N_edge, diffN, diffN_face,
316 diffN_edge, t_loc, NULL, NULL, dofs_x, dofs_x_edge, dofs_x_face,
317 &*kExtNodeFace.data().begin(), Kext_edge_face,
318 &*kExtFaceFace.data().begin(), gaussPts.size2(), &gaussPts(2, 0));
319 // cerr << "kExtNodeFace " << kExtNodeFace << endl;
320 // cerr << "kExtFaceFace " << kExtFaceFace << endl;
321 CHKERR MatSetValues(snes_B, 9, dofs_x_indices, kExtNodeFace.size2(),
322 dofs_x_face_indices, &*kExtNodeFace.data().begin(),
323 ADD_VALUES);
324 CHKERR MatSetValues(snes_B, kExtFaceFace.size1(), dofs_x_face_indices,
325 kExtFaceFace.size2(), dofs_x_face_indices,
326 &*kExtFaceFace.data().begin(), ADD_VALUES);
327 for (int ee = 0; ee < 3; ee++) {
328 // cerr << "kExtEdgeFace " << kExtEdgeFace[ee] << endl;
329 CHKERR MatSetValues(snes_B, kExtEdgeFace[ee].size1(),
330 dofs_x_edge_indices[ee], kExtFaceFace.size2(),
331 dofs_x_face_indices, Kext_edge_face[ee], ADD_VALUES);
332 }
333
334 kExtFaceEdge.resize(3);
335 kExtNodeEdge.resize(3);
336 kExtEdgeEdge.resize(3, 3);
337 for (int ee = 0; ee < 3; ee++) {
338 if (dOfs_x_edge_indices[ee].size() !=
339 (unsigned int)(3 * NBEDGE_H1(order_edge[ee]))) {
340 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
341 }
342 kExtFaceEdge[ee].resize(dOfs_x_face_indices.size(),
343 dOfs_x_edge_indices[ee].size());
344 kExtNodeEdge[ee].resize(9, dOfs_x_edge_indices[ee].size());
345 Kext_node_edge[ee] = &*kExtNodeEdge[ee].data().begin();
346 Kext_face_edge[ee] = &*kExtFaceEdge[ee].data().begin();
347 for (int EE = 0; EE < 3; EE++) {
348 kExtEdgeEdge(EE, ee).resize(dOfs_x_edge_indices[EE].size(),
349 dOfs_x_edge_indices[ee].size());
350 Kext_edge_edge[EE][ee] = &*kExtEdgeEdge(EE, ee).data().begin();
351 }
352 }
353 CHKERR KExt_hh_hierarchical_edge(
354 r * eps, order_face, order_edge, N, N_face, N_edge, diffN, diffN_face,
355 diffN_edge, t_loc, NULL, NULL, dofs_x, dofs_x_edge, dofs_x_face,
356 Kext_node_edge, Kext_edge_edge, Kext_face_edge, gaussPts.size2(),
357 &gaussPts(2, 0));
358 for (int ee = 0; ee < 3; ee++) {
359 CHKERR MatSetValues(snes_B, kExtFaceEdge[ee].size1(), dofs_x_face_indices,
360 kExtFaceEdge[ee].size2(), dofs_x_edge_indices[ee],
361 &*kExtFaceEdge[ee].data().begin(), ADD_VALUES);
362 CHKERR MatSetValues(snes_B, 9, dofs_x_indices, kExtNodeEdge[ee].size2(),
363 dofs_x_edge_indices[ee],
364 &*kExtNodeEdge[ee].data().begin(), ADD_VALUES);
365 for (int EE = 0; EE < 3; EE++) {
366 CHKERR MatSetValues(snes_B, kExtEdgeEdge(EE, ee).size1(),
367 dofs_x_edge_indices[EE], kExtEdgeEdge(EE, ee).size2(),
368 dofs_x_edge_indices[ee], Kext_edge_edge[EE][ee],
369 ADD_VALUES);
370 }
371 }
372
373 MoFEMFunctionReturn(0);
374}
375
376MoFEMErrorCode NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::
377 reBaseToFaceLoocalCoordSystem(MatrixDouble &t_glob_nodal) {
378 MoFEMFunctionBeginHot;
379 double s1[3], s2[3], normal[3], q[9];
380 CHKERR ShapeFaceBaseMBTRI(diffN, &*coords.data().begin(), normal, s1, s2);
381 double nrm2_normal = cblas_dnrm2(3, normal, 1);
382 cblas_dscal(3, 1. / nrm2_normal, normal, 1);
383 cblas_dcopy(3, s1, 1, &q[0], 1);
384 cblas_dcopy(3, s2, 1, &q[3], 1);
385 cblas_dcopy(3, normal, 1, &q[6], 1);
386 __CLPK_integer info;
387 __CLPK_integer ipiv[3];
388 info = lapack_dgesv(3, 3, q, 3, ipiv, &*t_glob_nodal.data().begin(), 3);
389 if (info != 0) {
390 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
391 "error solve dgesv info = %d", info);
392 }
393 MoFEMFunctionReturnHot(0);
394}
395
396MoFEMErrorCode
397NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::calcTraction() {
398 MoFEMFunctionBeginHot;
399 EntityHandle ent = numeredEntFiniteElementPtr->getEnt();
400 map<int, bCPressure>::iterator mip = mapPressure.begin();
401 tLoc.resize(3);
402 tLoc[0] = tLoc[1] = tLoc[2] = 0;
403 for (; mip != mapPressure.end(); mip++) {
404 if (mip->second.tRis.find(ent) != mip->second.tRis.end()) {
405 tLoc[2] -= mip->second.data.data.value1;
406 }
407 }
408 tLocNodal.resize(3, 3);
409 for (int nn = 0; nn < 3; nn++) {
410 for (int dd = 0; dd < 3; dd++) {
411 tLocNodal(nn, dd) = tLoc[dd];
412 }
413 }
414
415 map<int, bCForce>::iterator mif = mapForce.begin();
416 for (; mif != mapForce.end(); mif++) {
417 if (mif->second.tRis.find(ent) != mif->second.tRis.end()) {
418 tGlob.resize(3);
419 tGlob[0] = mif->second.data.data.value3;
420 tGlob[1] = mif->second.data.data.value4;
421 tGlob[2] = mif->second.data.data.value5;
422 tGlob *= mif->second.data.data.value1;
423 tGlobNodal.resize(3, 3);
424 for (int nn = 0; nn < 3; nn++) {
425 for (int dd = 0; dd < 3; dd++) {
426 tGlobNodal(nn, dd) = tGlob[dd];
427 }
428 }
429 CHKERR reBaseToFaceLoocalCoordSystem(tGlobNodal);
430 tLocNodal += tGlobNodal;
431 }
432 }
433
434 VectorDouble scale(1, 1);
435 CHKERR MethodForForceScaling::applyScale(this, methodsOp, scale);
436 tLocNodal *= scale[0];
437
438 t_loc = &*tLocNodal.data().begin();
439
440 MoFEMFunctionReturnHot(0);
441}
442
443MoFEMErrorCode
444NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::preProcess() {
445 MoFEMFunctionBeginHot;
446 PetscOptionsBegin(mField.get_comm(), "",
447 "Surface Pressure (complex for lazy)", "none");
448 PetscBool is_conservative = PETSC_TRUE;
449 CHKERR PetscOptionsBool("-is_conservative_force", "is conservative force", "",
450 PETSC_TRUE, &is_conservative, PETSC_NULLPTR);
451 if (is_conservative == PETSC_FALSE) {
452 typeOfForces = NONCONSERVATIVE;
453 }
454 PetscOptionsEnd();
455 MoFEMFunctionReturnHot(0);
456}
457
458MoFEMErrorCode NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::
459operator()() {
460 MoFEMFunctionBeginHot;
461
462 {
463
464 {
465
466 dofs_X = &*coords.data().begin();
467 for (int ee = 0; ee < 3; ee++) {
468 dofs_X_edge[ee] = NULL;
469 idofs_X_edge[ee] = NULL;
470 order_edge_material[ee] = 0;
471 }
472 dofs_X_face = NULL;
473 idofs_X_face = NULL;
474 order_face_material = 0;
475
476 dofs_x = &*coords.data().begin();
477 idofs_x = NULL;
478 for (int ee = 0; ee < 3; ee++) {
479 order_edge[ee] = 0;
480 N_edge[ee] = NULL;
481 diffN_edge[ee] = NULL;
482 dofs_x_edge[ee] = NULL;
483 idofs_x_edge[ee] = NULL;
484 }
485 order_face = 0;
486 N_face = NULL;
487 diffN_face = NULL;
488 dofs_x_face = NULL;
489 idofs_x_face = NULL;
490 }
491
492 CHKERR FaceElementForcesAndSourcesCore::operator()();
493 CHKERR calcTraction();
494
495 switch (snes_ctx) {
496 case CTX_SNESNONE:
497 case CTX_SNESSETFUNCTION: {
498 tLocNodal *= *sCaleRhs;
499 // cerr << "sCaleRhs " << *sCaleRhs << endl;
500 // cerr << tLocNodal << endl;
501 CHKERR rHs();
502 } break;
503 case CTX_SNESSETJACOBIAN: {
504 tLocNodal *= *sCaleLhs;
505 CHKERR lHs();
506 } break;
507 }
508 }
509
510 MoFEMFunctionReturnHot(0);
511}
512
513MoFEMErrorCode
514NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::addForce(int ms_id) {
515 MeshsetsManager *mesh_manager_ptr;
516 const CubitMeshSets *cubit_meshset_ptr;
517
518 MoFEMFunctionBeginHot;
519 CHKERR mField.getInterface(mesh_manager_ptr);
520 CHKERR mesh_manager_ptr->getCubitMeshsetPtr(ms_id, NODESET,
521 &cubit_meshset_ptr);
522 CHKERR cubit_meshset_ptr->getBcDataStructure(mapForce[ms_id].data);
523 CHKERR mField.get_moab().get_entities_by_type(
524 cubit_meshset_ptr->meshset, MBTRI, mapForce[ms_id].tRis, true);
525 MoFEMFunctionReturnHot(0);
526}
527
528MoFEMErrorCode
529NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::addPressure(
530 int ms_id) {
531 MeshsetsManager *mesh_manager_ptr;
532 const CubitMeshSets *cubit_meshset_ptr;
533
534 MoFEMFunctionBeginHot;
535 CHKERR mField.getInterface(mesh_manager_ptr);
536 CHKERR mesh_manager_ptr->getCubitMeshsetPtr(ms_id, SIDESET,
537 &cubit_meshset_ptr);
538 CHKERR cubit_meshset_ptr->getBcDataStructure(mapPressure[ms_id].data);
539 CHKERR mField.get_moab().get_entities_by_type(
540 cubit_meshset_ptr->meshset, MBTRI, mapPressure[ms_id].tRis, true);
541 MoFEMFunctionReturnHot(0);
542}
543
544NeumannForcesSurfaceComplexForLazy::NeumannForcesSurfaceComplexForLazy(
545 MoFEM::Interface &m_field, Mat _Aij, Vec _F, std::string spatial_field_name,
546 std::string material_field_name)
547 : mField(m_field), feSpatial(m_field, _Aij, _F, NULL, NULL,
548 spatial_field_name, material_field_name),
549 spatialField(spatial_field_name), materialField(material_field_name) {
550
551 double def_scale = 1.;
552 const EntityHandle root_meshset = mField.get_moab().get_root_set();
553 CHKERR mField.get_moab().tag_get_handle(
554 "_LoadFactor_Scale_", 1, MB_TYPE_DOUBLE, thScale,
555 MB_TAG_CREAT | MB_TAG_EXCL | MB_TAG_MESH, &def_scale);
556 if (rval == MB_ALREADY_ALLOCATED) {
557 rval = mField.get_moab().tag_get_by_ptr(thScale, &root_meshset, 1,
558 (const void **)&sCale);
559 MOAB_THROW(rval);
560 } else {
561 MOAB_THROW(rval);
562 rval =
563 mField.get_moab().tag_set_data(thScale, &root_meshset, 1, &def_scale);
564 MOAB_THROW(rval);
565 rval = mField.get_moab().tag_get_by_ptr(thScale, &root_meshset, 1,
566 (const void **)&sCale);
567 MOAB_THROW(rval);
568 }
569
570 feSpatial.sCaleLhs = sCale;
571 feSpatial.sCaleRhs = sCale;
572}
573
574NeumannForcesSurfaceComplexForLazy::NeumannForcesSurfaceComplexForLazy(
575 MoFEM::Interface &m_field, Mat _Aij, Vec _F, double *scale_lhs,
576 double *scale_rhs, std::string spatial_field_name,
577 std::string material_field_name)
578 : mField(m_field), feSpatial(m_field, _Aij, _F, scale_lhs, scale_rhs,
579 spatial_field_name, material_field_name),
580 spatialField(spatial_field_name), materialField(material_field_name) {}
581
582extern "C" {
583
584// External forces
585MoFEMErrorCode Traction_hierarchical(int order, int *order_edge, double *N,
586 double *N_face, double *N_edge[],
587 double *t, double *t_edge[],
588 double *t_face, double *traction, int gg) {
589 MoFEMFunctionBeginHot;
590 int dd, ee;
591 for (dd = 0; dd < 3; dd++)
592 traction[dd] = cblas_ddot(3, &N[gg * 3], 1, &t[dd], 3);
593 if (t_face != NULL) {
594 int nb_dofs_face = NBFACETRI_H1(order);
595 if (nb_dofs_face > 0) {
596 for (dd = 0; dd < 3; dd++)
597 traction[dd] += cblas_ddot(nb_dofs_face, &N_face[gg * nb_dofs_face], 1,
598 &t_face[dd], 3);
599 }
600 }
601 if (t_edge != NULL) {
602 ee = 0;
603 for (; ee < 3; ee++) {
604 if (t_edge[ee] == NULL)
605 continue;
606 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
607 if (nb_dofs_edge > 0) {
608 for (dd = 0; dd < 3; dd++) {
609 traction[dd] +=
610 cblas_ddot(nb_dofs_edge, &(N_edge[ee][gg * nb_dofs_edge]), 1,
611 &(t_edge[ee][dd]), 3);
612 }
613 }
614 }
615 }
616 MoFEMFunctionReturnHot(0);
617}
618MoFEMErrorCode Fext_h_hierarchical(
619 int order, int *order_edge, double *N, double *N_face, double *N_edge[],
620 double *diffN, double *diffN_face, double *diffN_edge[], double *t,
621 double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[],
622 double *dofs_x_face, double *idofs_x, double *idofs_x_edge[],
623 double *idofs_x_face, double *Fext, double *Fext_edge[], double *Fext_face,
624 double *iFext, double *iFext_edge[], double *iFext_face, int g_dim,
625 const double *g_w) {
626 MoFEMFunctionBeginHot;
627 int dd, nn, ee, gg;
628 if (Fext != NULL)
629 bzero(Fext, 9 * sizeof(double));
630 if (iFext != NULL)
631 bzero(iFext, 9 * sizeof(double));
632 ee = 0;
633 for (; ee < 3; ee++) {
634 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
635 if (nb_dofs_edge == 0)
636 continue;
637 if (Fext_edge != NULL)
638 bzero(Fext_edge[ee], 3 * nb_dofs_edge * sizeof(double));
639 if (iFext_edge != NULL)
640 bzero(iFext_edge[ee], 3 * nb_dofs_edge * sizeof(double));
641 }
642 int nb_dofs_face = NBFACETRI_H1(order);
643 if (nb_dofs_face != 0) {
644 if (Fext_face != NULL)
645 bzero(Fext_face, 3 * nb_dofs_face * sizeof(double));
646 if (iFext_face != NULL)
647 bzero(iFext_face, 3 * nb_dofs_face * sizeof(double));
648 }
649 gg = 0;
650 for (; gg < g_dim; gg++) {
651 double traction[3] = {0, 0, 0};
652 CHKERR Traction_hierarchical(order, order_edge, N, N_face, N_edge, t,
653 t_edge, t_face, traction, gg);
654 __CLPK_doublecomplex xnormal[3], xs1[3], xs2[3];
655 CHKERR Normal_hierarchical(
656 // FIXME: order of tractions approximation could be different for field
657 // approx.
658 order, order_edge, order, order_edge, diffN, diffN_face, diffN_edge,
659 dofs_x, dofs_x_edge, dofs_x_face, idofs_x, idofs_x_edge, idofs_x_face,
660 xnormal, xs1, xs2, gg);
661 CHKERR Base_scale(xnormal, xs1, xs2);
662 double normal_real[3], s1_real[3], s2_real[3];
663 double normal_imag[3], s1_imag[3], s2_imag[3];
664 for (dd = 0; dd < 3; dd++) {
665 normal_real[dd] = 0.5 * xnormal[dd].r;
666 normal_imag[dd] = 0.5 * xnormal[dd].i;
667 s1_real[dd] = 0.5 * xs1[dd].r;
668 s1_imag[dd] = 0.5 * xs1[dd].i;
669 s2_real[dd] = 0.5 * xs2[dd].r;
670 s2_imag[dd] = 0.5 * xs2[dd].i;
671 }
672 nn = 0;
673 for (; nn < 3; nn++) {
674 if (Fext != NULL)
675 for (dd = 0; dd < 3; dd++) {
676 // fprintf(stderr,"%d %f %f %f %f %f %f\n",
677 // gg,g_w[gg],N[3*gg+nn],normal_real[dd],
678 // traction[0],traction[1],traction[2]);
679 Fext[3 * nn + dd] +=
680 g_w[gg] * N[3 * gg + nn] * normal_real[dd] * traction[2];
681 Fext[3 * nn + dd] +=
682 g_w[gg] * N[3 * gg + nn] * s1_real[dd] * traction[0];
683 Fext[3 * nn + dd] +=
684 g_w[gg] * N[3 * gg + nn] * s2_real[dd] * traction[1];
685 }
686 if (iFext != NULL)
687 for (dd = 0; dd < 3; dd++) {
688 iFext[3 * nn + dd] +=
689 g_w[gg] * N[3 * gg + nn] * normal_imag[dd] * traction[2];
690 iFext[3 * nn + dd] +=
691 g_w[gg] * N[3 * gg + nn] * s1_imag[dd] * traction[0];
692 iFext[3 * nn + dd] +=
693 g_w[gg] * N[3 * gg + nn] * s2_imag[dd] * traction[1];
694 }
695 }
696 ee = 0;
697 for (; ee < 3; ee++) {
698 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
699 if (nb_dofs_edge == 0)
700 continue;
701 int nn = 0;
702 for (; nn < nb_dofs_edge; nn++) {
703 if (Fext_edge != NULL)
704 for (dd = 0; dd < 3; dd++) {
705 Fext_edge[ee][3 * nn + dd] += g_w[gg] *
706 N_edge[ee][gg * nb_dofs_edge + nn] *
707 normal_real[dd] * traction[2];
708 Fext_edge[ee][3 * nn + dd] += g_w[gg] *
709 N_edge[ee][gg * nb_dofs_edge + nn] *
710 s1_real[dd] * traction[0];
711 Fext_edge[ee][3 * nn + dd] += g_w[gg] *
712 N_edge[ee][gg * nb_dofs_edge + nn] *
713 s2_real[dd] * traction[1];
714 }
715 if (iFext_edge != NULL) {
716 for (dd = 0; dd < 3; dd++) {
717 iFext_edge[ee][3 * nn + dd] += g_w[gg] *
718 N_edge[ee][gg * nb_dofs_edge + nn] *
719 normal_imag[dd] * traction[2];
720 iFext_edge[ee][3 * nn + dd] += g_w[gg] *
721 N_edge[ee][gg * nb_dofs_edge + nn] *
722 s1_imag[dd] * traction[0];
723 iFext_edge[ee][3 * nn + dd] += g_w[gg] *
724 N_edge[ee][gg * nb_dofs_edge + nn] *
725 s2_imag[dd] * traction[1];
726 }
727 }
728 }
729 }
730 if (nb_dofs_face != 0) {
731 nn = 0;
732 for (; nn < nb_dofs_face; nn++) {
733 if (Fext_face != NULL)
734 for (dd = 0; dd < 3; dd++) {
735 Fext_face[3 * nn + dd] += g_w[gg] * N_face[gg * nb_dofs_face + nn] *
736 normal_real[dd] * traction[2];
737 Fext_face[3 * nn + dd] += g_w[gg] * N_face[gg * nb_dofs_face + nn] *
738 s1_real[dd] * traction[0];
739 Fext_face[3 * nn + dd] += g_w[gg] * N_face[gg * nb_dofs_face + nn] *
740 s2_real[dd] * traction[1];
741 }
742 if (iFext_face != NULL)
743 for (dd = 0; dd < 3; dd++) {
744 iFext_face[3 * nn + dd] += g_w[gg] *
745 N_face[gg * nb_dofs_face + nn] *
746 normal_imag[dd] * traction[2];
747 iFext_face[3 * nn + dd] += g_w[gg] *
748 N_face[gg * nb_dofs_face + nn] *
749 s1_imag[dd] * traction[0];
750 iFext_face[3 * nn + dd] += g_w[gg] *
751 N_face[gg * nb_dofs_face + nn] *
752 s2_imag[dd] * traction[1];
753 }
754 }
755 }
756 }
757 MoFEMFunctionReturnHot(0);
758}
759MoFEMErrorCode KExt_hh_hierarchical(
760 double eps, int order, int *order_edge, double *N, double *N_face,
761 double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[],
762 double *t, double *t_edge[], double *t_face, double *dofs_x,
763 double *dofs_x_edge[], double *dofs_x_face, double *KExt_hh,
764 double *KExt_edgeh[], double *KExt_faceh, int g_dim, const double *g_w) {
765 MoFEMFunctionBeginHot;
766 int gg, dd, ii, nn, ee;
767 bzero(KExt_hh, 9 * 9 * sizeof(double));
768 if (KExt_edgeh != NULL) {
769 for (ee = 0; ee < 3; ee++) {
770 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
771 bzero(KExt_edgeh[ee], 9 * 3 * nb_dofs_edge * sizeof(double));
772 }
773 }
774 int nb_dofs_face = NBFACETRI_H1(order);
775 if (KExt_faceh != NULL) {
776 bzero(KExt_faceh, 9 * 3 * nb_dofs_face * sizeof(double));
777 }
778 for (gg = 0; gg < g_dim; gg++) {
779 double traction[3] = {0, 0, 0};
780 CHKERR Traction_hierarchical(order, order_edge, N, N_face, N_edge, t,
781 t_edge, t_face, traction, gg);
782 //
783 __CLPK_doublecomplex xnormal[3], xs1[3], xs2[3];
784 double idofs_x[9];
785 for (ii = 0; ii < 9; ii++) {
786 bzero(idofs_x, 9 * sizeof(double));
787 idofs_x[ii] = eps;
788 CHKERR Normal_hierarchical(order, order_edge, // FIXME
789 order, order_edge, diffN, diffN_face,
790 diffN_edge, dofs_x, dofs_x_edge, dofs_x_face,
791 idofs_x, NULL, NULL, xnormal, xs1, xs2, gg);
792 CHKERR Base_scale(xnormal, xs1, xs2);
793 double normal_imag[3], s1_imag[3], s2_imag[3];
794 for (dd = 0; dd < 3; dd++) {
795 normal_imag[dd] = 0.5 * xnormal[dd].i / eps;
796 s1_imag[dd] = 0.5 * xs1[dd].i / eps;
797 s2_imag[dd] = 0.5 * xs2[dd].i / eps;
798 }
799 nn = 0;
800 for (; nn < 3; nn++) {
801 for (dd = 0; dd < 3; dd++) {
802 KExt_hh[ii + 9 * 3 * nn + 9 * dd] +=
803 g_w[gg] * N[3 * gg + nn] * normal_imag[dd] * traction[2];
804 KExt_hh[ii + 9 * 3 * nn + 9 * dd] +=
805 g_w[gg] * N[3 * gg + nn] * s1_imag[dd] * traction[0];
806 KExt_hh[ii + 9 * 3 * nn + 9 * dd] +=
807 g_w[gg] * N[3 * gg + nn] * s2_imag[dd] * traction[1];
808 }
809 }
810 if (KExt_edgeh != NULL) {
811 for (ee = 0; ee < 3; ee++) {
812 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
813 for (nn = 0; nn < nb_dofs_edge; nn++) {
814 for (dd = 0; dd < 3; dd++) {
815 KExt_edgeh[ee][ii + 9 * 3 * nn + 9 * dd] +=
816 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] *
817 normal_imag[dd] * traction[2];
818 KExt_edgeh[ee][ii + 9 * 3 * nn + 9 * dd] +=
819 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s1_imag[dd] *
820 traction[0];
821 KExt_edgeh[ee][ii + 9 * 3 * nn + 9 * dd] +=
822 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s2_imag[dd] *
823 traction[1];
824 }
825 }
826 }
827 }
828 if (KExt_faceh != NULL) {
829 for (nn = 0; nn < nb_dofs_face; nn++) {
830 for (dd = 0; dd < 3; dd++) {
831 KExt_faceh[ii + 3 * 9 * nn + 9 * dd] +=
832 g_w[gg] * N_face[nb_dofs_face * gg + nn] * normal_imag[dd] *
833 traction[2];
834 KExt_faceh[ii + 3 * 9 * nn + 9 * dd] +=
835 g_w[gg] * N_face[nb_dofs_face * gg + nn] * s1_imag[dd] *
836 traction[0];
837 KExt_faceh[ii + 3 * 9 * nn + 9 * dd] +=
838 g_w[gg] * N_face[nb_dofs_face * gg + nn] * s2_imag[dd] *
839 traction[1];
840 }
841 }
842 }
843 }
844 }
845 MoFEMFunctionReturnHot(0);
846}
847MoFEMErrorCode KExt_hh_hierarchical_edge(
848 double eps, int order, int *order_edge, double *N, double *N_face,
849 double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[],
850 double *t, double *t_edge[], double *t_face, double *dofs_x,
851 double *dofs_x_edge[], double *dofs_x_face, double *KExt_hedge[3],
852 double *KExt_edgeedge[3][3], double *KExt_faceedge[3], int g_dim,
853 const double *g_w) {
854 MoFEMFunctionBeginHot;
855 int gg, dd, ii, nn, ee, EE;
856 int nb_dofs_face = NBFACETRI_H1(order);
857 for (EE = 0; EE < 3; EE++) {
858 int nb_dofs_edge_EE = NBEDGE_H1(order_edge[EE]);
859 bzero(KExt_hedge[EE], 9 * 3 * nb_dofs_edge_EE * sizeof(double));
860 if (KExt_edgeedge != NULL) {
861 for (ee = 0; ee < 3; ee++) {
862 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
863 bzero(KExt_edgeedge[EE][ee],
864 3 * nb_dofs_edge_EE * 3 * nb_dofs_edge * sizeof(double));
865 }
866 }
867 if (KExt_faceedge != NULL) {
868 bzero(KExt_faceedge[EE],
869 3 * nb_dofs_edge_EE * 3 * nb_dofs_face * sizeof(double));
870 }
871 }
872 for (gg = 0; gg < g_dim; gg++) {
873 double traction[3] = {0, 0, 0};
874 CHKERR Traction_hierarchical(order, order_edge, N, N_face, N_edge, t,
875 t_edge, t_face, traction, gg);
876 for (EE = 0; EE < 3; EE++) {
877 int nb_dofs_edge_EE = NBEDGE_H1(order_edge[EE]);
878 double *idofs_x_edge[3] = {NULL, NULL, NULL};
879 double idofs_x_edge_EE[3 * nb_dofs_edge_EE];
880 idofs_x_edge[EE] = idofs_x_edge_EE;
881 for (ii = 0; ii < 3 * nb_dofs_edge_EE; ii++) {
882 bzero(idofs_x_edge_EE, 3 * nb_dofs_edge_EE * sizeof(double));
883 idofs_x_edge_EE[ii] = eps;
884 __CLPK_doublecomplex xnormal[3], xs1[3], xs2[3];
885 CHKERR Normal_hierarchical(order, order_edge, // FIXME
886 order, order_edge, diffN, diffN_face,
887 diffN_edge, dofs_x, dofs_x_edge, dofs_x_face,
888 NULL, idofs_x_edge, NULL, xnormal, xs1, xs2,
889 gg);
890 CHKERR Base_scale(xnormal, xs1, xs2);
891 double normal_imag[3], s1_imag[3], s2_imag[3];
892 for (dd = 0; dd < 3; dd++) {
893 normal_imag[dd] = 0.5 * xnormal[dd].i / eps;
894 s1_imag[dd] = 0.5 * xs1[dd].i / eps;
895 s2_imag[dd] = 0.5 * xs2[dd].i / eps;
896 }
897 for (nn = 0; nn < 3; nn++) {
898 for (dd = 0; dd < 3; dd++) {
899 KExt_hedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
900 3 * nb_dofs_edge_EE * dd] +=
901 g_w[gg] * N[3 * gg + nn] * normal_imag[dd] * traction[2];
902 KExt_hedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
903 3 * nb_dofs_edge_EE * dd] +=
904 g_w[gg] * N[3 * gg + nn] * s1_imag[dd] * traction[0];
905 KExt_hedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
906 3 * nb_dofs_edge_EE * dd] +=
907 g_w[gg] * N[3 * gg + nn] * s2_imag[dd] * traction[1];
908 }
909 }
910 if (KExt_edgeedge != NULL) {
911 for (ee = 0; ee < 3; ee++) {
912 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
913 for (nn = 0; nn < nb_dofs_edge; nn++) {
914 for (dd = 0; dd < 3; dd++) {
915 KExt_edgeedge[EE][ee][ii + 3 * nb_dofs_edge_EE * 3 * nn +
916 3 * nb_dofs_edge_EE * dd] +=
917 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] *
918 normal_imag[dd] * traction[2];
919 KExt_edgeedge[EE][ee][ii + 3 * nb_dofs_edge_EE * 3 * nn +
920 3 * nb_dofs_edge_EE * dd] +=
921 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s1_imag[dd] *
922 traction[0];
923 KExt_edgeedge[EE][ee][ii + 3 * nb_dofs_edge_EE * 3 * nn +
924 3 * nb_dofs_edge_EE * dd] +=
925 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s2_imag[dd] *
926 traction[1];
927 }
928 }
929 }
930 }
931 if (KExt_faceedge != NULL) {
932 for (nn = 0; nn < nb_dofs_face; nn++) {
933 for (dd = 0; dd < 3; dd++) {
934 KExt_faceedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
935 3 * nb_dofs_edge_EE * dd] +=
936 g_w[gg] * N_face[nb_dofs_face * gg + nn] * normal_imag[dd] *
937 traction[2];
938 KExt_faceedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
939 3 * nb_dofs_edge_EE * dd] +=
940 g_w[gg] * N_face[nb_dofs_face * gg + nn] * s1_imag[dd] *
941 traction[0];
942 KExt_faceedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
943 3 * nb_dofs_edge_EE * dd] +=
944 g_w[gg] * N_face[nb_dofs_face * gg + nn] * s2_imag[dd] *
945 traction[1];
946 }
947 }
948 }
949 }
950 }
951 }
952 MoFEMFunctionReturnHot(0);
953}
954MoFEMErrorCode
955KExt_hh_hierarchical_face(double eps, int order, int *order_edge, double *N,
956 double *N_face, double *N_edge[], double *diffN,
957 double *diffN_face, double *diffN_edge[], double *t,
958 double *t_edge[], double *t_face, double *dofs_x,
959 double *dofs_x_edge[], double *dofs_x_face,
960 double *KExt_hface, double *KExt_edgeface[3],
961 double *KExt_faceface, int g_dim, const double *g_w) {
962 MoFEMFunctionBeginHot;
963 int gg, dd, ii, nn, ee;
964 int nb_dofs_face = NBFACETRI_H1(order);
965 bzero(KExt_hface, 9 * 3 * nb_dofs_face * sizeof(double));
966 if (KExt_edgeface != NULL) {
967 for (ee = 0; ee < 3; ee++) {
968 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
969 bzero(KExt_edgeface[ee],
970 3 * nb_dofs_face * 3 * nb_dofs_edge * sizeof(double));
971 }
972 }
973 if (KExt_faceface != NULL) {
974 bzero(KExt_faceface, 3 * nb_dofs_face * 3 * nb_dofs_face * sizeof(double));
975 }
976 for (gg = 0; gg < g_dim; gg++) {
977 double traction[3] = {0, 0, 0};
978 CHKERR Traction_hierarchical(order, order_edge, N, N_face, N_edge, t,
979 t_edge, t_face, traction, gg);
980 double idofs_x_face[3 * nb_dofs_face];
981 for (ii = 0; ii < 3 * nb_dofs_face; ii++) {
982 bzero(idofs_x_face, 3 * nb_dofs_face * sizeof(double));
983 idofs_x_face[ii] = eps;
984 __CLPK_doublecomplex xnormal[3], xs1[3], xs2[3];
985 CHKERR Normal_hierarchical(
986 order, order_edge, // FIXME
987 order, order_edge, diffN, diffN_face, diffN_edge, dofs_x, dofs_x_edge,
988 dofs_x_face, NULL, NULL, idofs_x_face, xnormal, xs1, xs2, gg);
989 CHKERR Base_scale(xnormal, xs1, xs2);
990 double normal_imag[3], s1_imag[3], s2_imag[3];
991 for (dd = 0; dd < 3; dd++) {
992 normal_imag[dd] = 0.5 * xnormal[dd].i / eps;
993 s1_imag[dd] = 0.5 * xs1[dd].i / eps;
994 s2_imag[dd] = 0.5 * xs2[dd].i / eps;
995 }
996 for (nn = 0; nn < 3; nn++) {
997 for (dd = 0; dd < 3; dd++) {
998 KExt_hface[ii + 3 * nb_dofs_face * 3 * nn + 3 * nb_dofs_face * dd] +=
999 g_w[gg] * N[3 * gg + nn] * normal_imag[dd] * traction[2];
1000 KExt_hface[ii + 3 * nb_dofs_face * 3 * nn + 3 * nb_dofs_face * dd] +=
1001 g_w[gg] * N[3 * gg + nn] * s1_imag[dd] * traction[0];
1002 KExt_hface[ii + 3 * nb_dofs_face * 3 * nn + 3 * nb_dofs_face * dd] +=
1003 g_w[gg] * N[3 * gg + nn] * s2_imag[dd] * traction[1];
1004 }
1005 }
1006 if (KExt_edgeface != NULL) {
1007 for (ee = 0; ee < 3; ee++) {
1008 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
1009 for (nn = 0; nn < nb_dofs_edge; nn++) {
1010 for (dd = 0; dd < 3; dd++) {
1011 KExt_edgeface[ee][ii + 3 * nb_dofs_face * 3 * nn +
1012 3 * nb_dofs_face * dd] +=
1013 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] *
1014 normal_imag[dd] * traction[2];
1015 KExt_edgeface[ee][ii + 3 * nb_dofs_face * 3 * nn +
1016 3 * nb_dofs_face * dd] +=
1017 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s1_imag[dd] *
1018 traction[0];
1019 KExt_edgeface[ee][ii + 3 * nb_dofs_face * 3 * nn +
1020 3 * nb_dofs_face * dd] +=
1021 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s2_imag[dd] *
1022 traction[1];
1023 }
1024 }
1025 }
1026 }
1027 if (KExt_faceface != NULL) {
1028 for (nn = 0; nn < nb_dofs_face; nn++) {
1029 for (dd = 0; dd < 3; dd++) {
1030 KExt_faceface[ii + 3 * nb_dofs_face * 3 * nn +
1031 3 * nb_dofs_face * dd] +=
1032 g_w[gg] * N_face[nb_dofs_face * gg + nn] * normal_imag[dd] *
1033 traction[2];
1034 KExt_faceface[ii + 3 * nb_dofs_face * 3 * nn +
1035 3 * nb_dofs_face * dd] +=
1036 g_w[gg] * N_face[nb_dofs_face * gg + nn] * s1_imag[dd] *
1037 traction[0];
1038 KExt_faceface[ii + 3 * nb_dofs_face * 3 * nn +
1039 3 * nb_dofs_face * dd] +=
1040 g_w[gg] * N_face[nb_dofs_face * gg + nn] * s2_imag[dd] *
1041 traction[1];
1042 }
1043 }
1044 }
1045 }
1046 }
1047 MoFEMFunctionReturnHot(0);
1048}
1049}
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition HookeElement.hpp:75
Fext_h_hierarchical
MoFEMErrorCode Fext_h_hierarchical(int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[], double *t, double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[], double *dofs_x_face, double *idofs_x, double *idofs_x_edge[], double *idofs_x_face, double *Fext, double *Fext_egde[], double *Fext_face, double *iFext, double *iFext_egde[], double *iFext_face, int g_dim, const double *g_w)
Definition SurfacePressureComplexForLazy.cpp:618
KExt_hh_hierarchical
MoFEMErrorCode KExt_hh_hierarchical(double eps, int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[], double *t, double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[], double *dofs_x_face, double *KExt_hh, double *KExt_egdeh[3], double *KExt_faceh, int g_dim, const double *g_w)
tricircumcenter3d_tp
void tricircumcenter3d_tp(double a[3], double b[3], double c[3], double circumcenter[3], double *xi, double *eta)
Traction_hierarchical
MoFEMErrorCode Traction_hierarchical(int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *t, double *t_edge[], double *t_face, double *traction, int gg)
Definition SurfacePressureComplexForLazy.cpp:585
KExt_hh_hierarchical_edge
MoFEMErrorCode KExt_hh_hierarchical_edge(double eps, int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[], double *t, double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[], double *dofs_x_face, double *Khext_edge[3], double *KExt_edgeegde[3][3], double *KExt_faceedge[3], int g_dim, const double *g_w)
Definition SurfacePressureComplexForLazy.cpp:847
KExt_hh_hierarchical_face
MoFEMErrorCode KExt_hh_hierarchical_face(double eps, int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[], double *t, double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[], double *dofs_x_face, double *KExt_hface, double *KExt_egdeface[3], double *KExt_faceface, int g_dim, const double *g_w)
Definition SurfacePressureComplexForLazy.cpp:955
tetcircumcenter_tp
void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3], double circumcenter[3], double *xi, double *eta, double *zeta)
a
constexpr double a
Definition approx_sphere.cpp:30
eps
static const double eps
Definition check_base_functions_derivatives_on_tet.cpp:11
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition definitions.h:557
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
H1
@ H1
continuous field
Definition definitions.h:85
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
NODESET
@ NODESET
Definition definitions.h:159
SIDESET
@ SIDESET
Definition definitions.h:160
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
order
constexpr int order
Definition dg_projection.cpp:18
ShapeFaceBaseMBTRI
PetscErrorCode ShapeFaceBaseMBTRI(double *diffN, const double *coords, double *normal, double *s1, double *s2)
Definition fem_tools.c:204
Base_scale
PetscErrorCode Base_scale(__CLPK_doublecomplex *xnormal, __CLPK_doublecomplex *xs1, __CLPK_doublecomplex *xs2)
Definition fem_tools.c:741
Normal_hierarchical
PetscErrorCode Normal_hierarchical(int order_approx, int *order_edge_approx, int order, int *order_edge, double *diffN, double *diffN_face, double *diffN_edge[], double *dofs, double *dofs_edge[], double *dofs_face, double *idofs, double *idofs_edge[], double *idofs_face, __CLPK_doublecomplex *xnormal, __CLPK_doublecomplex *s1, __CLPK_doublecomplex *s2, int gg)
Complex normal.
Definition fem_tools.c:569
F
@ F
Definition free_surface.cpp:396
eta
double eta
Definition free_surface.cpp:172
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition MeshsetsManager.cpp:589
NBEDGE_H1
#define NBEDGE_H1(P)
Number of base function on edge for H1 space.
Definition h1_hdiv_hcurl_l2.h:55
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition h1_hdiv_hcurl_l2.h:60
c
const double c
speed of light (cm/ns)
Definition initial_diffusion.cpp:39
lapack_dgesv
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition lapack_wrap.h:188
__CLPK_integer
long int __CLPK_integer
Definition lapack_wrap.h:23
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition Exceptions.hpp:74
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition EntitiesFieldData.hpp:1666
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition EntitiesFieldData.hpp:1611
scale
double scale
Definition plastic.cpp:123
zeta
double zeta
Viscous hardening.
Definition plastic.cpp:130
t
constexpr double t
plate stiffness
Definition plate.cpp:58
field_name
constexpr auto field_name
Definition poisson_2d_homogeneous.cpp:13
_F
#define _F(n)
Definition quad.c:25
N
const int N
Definition speed_test.cpp:3
MethodForForceScaling::applyScale
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Definition MethodForForceScaling.hpp:21
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CubitMeshSets
this struct keeps basic methods for moab meshset about material and boundary conditions
Definition BCMultiIndices.hpp:19
MoFEM::CubitMeshSets::getBcDataStructure
MoFEMErrorCode getBcDataStructure(CUBIT_BC_DATA_TYPE &data) const
Definition BCMultiIndices.hpp:296
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition EntitiesFieldData.hpp:128
MoFEM::EntitiesFieldData::EntData::getOrder
ApproximationOrder getOrder() const
get approximation order
Definition EntitiesFieldData.hpp:1210
MoFEM::EntitiesFieldData::EntData::getDiffN
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
Definition EntitiesFieldData.hpp:1345
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition EntitiesFieldData.hpp:1318
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition EntitiesFieldData.hpp:1254
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition EntitiesFieldData.hpp:1214
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition ForcesAndSourcesCore.hpp:83
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition MeshsetsManager.hpp:104
MoFEM::SnesMethod::snes_f
Vec & snes_f
residual
Definition LoopMethods.hpp:127
MoFEM::SnesMethod::snes_B
Mat & snes_B
preconditioner of jacobian matrix
Definition LoopMethods.hpp:129
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
NeumannForcesSurfaceComplexForLazy::AuxMethodMaterial::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition SurfacePressureComplexForLazy.cpp:117
NeumannForcesSurfaceComplexForLazy::AuxMethodMaterial::AuxMethodMaterial
AuxMethodMaterial(const string &field_name, MyTriangleSpatialFE *my_ptr, const char type)
Definition SurfacePressureComplexForLazy.cpp:61
NeumannForcesSurfaceComplexForLazy::AuxMethodSpatial::AuxMethodSpatial
AuxMethodSpatial(const string &field_name, MyTriangleSpatialFE *my_ptr, const char type)
Definition SurfacePressureComplexForLazy.cpp:56
NeumannForcesSurfaceComplexForLazy::AuxMethodSpatial::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition SurfacePressureComplexForLazy.cpp:66
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::reBaseToFaceLoocalCoordSystem
MoFEMErrorCode reBaseToFaceLoocalCoordSystem(MatrixDouble &t_glob_nodal)
Definition SurfacePressureComplexForLazy.cpp:377
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition SurfacePressureComplexForLazy.cpp:459
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition SurfacePressureComplexForLazy.cpp:444
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::MyTriangleSpatialFE
MyTriangleSpatialFE(MoFEM::Interface &_mField, Mat _Aij, Vec &_F, double *scale_lhs, double *scale_rhs, std::string spatial_field_name="SPATIAL_POSITION", std::string mat_field_name="MESH_NODE_POSITIONS")
Definition SurfacePressureComplexForLazy.cpp:161
NeumannForcesSurfaceComplexForLazy::NeumannForcesSurfaceComplexForLazy
NeumannForcesSurfaceComplexForLazy(MoFEM::Interface &m_field, Mat _Aij, Vec _F, double *scale_lhs, double *scale_rhs, std::string spatial_field_name="SPATIAL_POSITION", std::string material_field_name="MESH_NODE_POSITIONS")
Definition SurfacePressureComplexForLazy.cpp:574
__CLPK_doublecomplex
Definition lapack_wrap.h:34
__CLPK_doublecomplex::i
__CLPK_doublereal i
Definition lapack_wrap.h:35
__CLPK_doublecomplex::r
__CLPK_doublereal r
Definition lapack_wrap.h:35
tricircumcenter3d_tp
void tricircumcenter3d_tp(a, b, c, circumcenter, double *xi, double *eta)
Definition tricircumcenter.c:100
Generated by Doxygen 1.12.0 and hosted at University of Glasgow