v0.14.0
Loading...
Searching...
No Matches
SurfacePressure.cpp
Go to the documentation of this file.
1/* \file SurfacePressure.cpp
2 \brief Implementation of pressure and forces on triangles surface
3
4 \todo Note that ALE version was not tested for quad elements, and is not
5 guarantee that will work.
6*/
7
8
9
10using namespace MoFEM;
11
12using namespace boost::numeric;
13
14MoFEMErrorCode NeumannForcesSurface::LinearVaringPresssure::getForce(
15 const EntityHandle ent, const VectorDouble3 &coords,
16 const VectorDouble3 &normal, VectorDouble3 &force) {
17 MoFEMFunctionBegin;
18 const double p = inner_prod(coords, linearConstants) + pressureShift;
19 force = normal * p / norm_2(normal);
20 MoFEMFunctionReturn(0);
21}
22
23NeumannForcesSurface::MyTriangleFE::MyTriangleFE(MoFEM::Interface &m_field)
24 : FaceElementForcesAndSourcesCore(m_field), addToRule(1) {}
25
26NeumannForcesSurface::OpNeumannForce::OpNeumannForce(
27 const std::string field_name, Vec _F, bCForce &data,
28 boost::ptr_vector<MethodForForceScaling> &methods_op, bool ho_geometry)
29 : FaceElementForcesAndSourcesCore::UserDataOperator(
30 field_name, UserDataOperator::OPROW),
31 F(_F), dAta(data), methodsOp(methods_op), hoGeometry(ho_geometry) {}
32
33MoFEMErrorCode
34NeumannForcesSurface::OpNeumannForce::doWork(int side, EntityType type,
35 EntitiesFieldData::EntData &data) {
36 MoFEMFunctionBegin;
37
38 if (data.getIndices().size() == 0)
39 MoFEMFunctionReturnHot(0);
40 EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
41 if (dAta.tRis.find(ent) == dAta.tRis.end())
42 MoFEMFunctionReturnHot(0);
43
44 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
45 int nb_row_dofs = data.getIndices().size() / rank;
46
47 Nf.resize(data.getIndices().size(), false);
48 Nf.clear();
49
50 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
51
52 for (unsigned int gg = 0; gg != data.getN().size1(); ++gg) {
53
54 // get integration weight and Jacobian of integration point (area of face)
55 double val = getGaussPts()(2, gg);
56 if (hoGeometry || fe_type == MBQUAD) {
57 val *= cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
58 if (fe_type == MBTRI)
59 val /= 2;
60 } else
61 val *= getArea();
62
63 // use data from module
64 for (int rr = 0; rr != rank; ++rr) {
65 double force;
66 if (rr == 0)
67 force = dAta.data.data.value3;
68 else if (rr == 1)
69 force = dAta.data.data.value4;
70 else if (rr == 2)
71 force = dAta.data.data.value5;
72 else
73 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
74 "data inconsistency");
75
76 force *= dAta.data.data.value1;
77 cblas_daxpy(nb_row_dofs, val * force, &data.getN()(gg, 0), 1, &Nf[rr],
78 rank);
79 }
80 }
81
82 // Scale force using user defined scaling operator
83 CHKERR MethodForForceScaling::applyScale(getFEMethod(), methodsOp, Nf);
84
85 auto get_f = [&]() {
86 if (F == PETSC_NULL)
87 return getKSPf();
88 return F;
89 };
90
91 // Assemble force into vector
92 CHKERR VecSetValues(get_f(), data, &*Nf.data().begin(), ADD_VALUES);
93
94 MoFEMFunctionReturn(0);
95}
96
97NeumannForcesSurface::OpNeumannForceAnalytical::OpNeumannForceAnalytical(
98 const std::string field_name, Vec f, const Range tris,
99 boost::ptr_vector<MethodForForceScaling> &methods_op,
100 boost::shared_ptr<MethodForAnalyticalForce> &analytical_force_op,
101 const bool ho_geometry)
102 : FaceElementForcesAndSourcesCore::UserDataOperator(
103 field_name, UserDataOperator::OPROW),
104 F(f), tRis(tris), methodsOp(methods_op),
105 analyticalForceOp(analytical_force_op), hoGeometry(ho_geometry) {}
106
107MoFEMErrorCode NeumannForcesSurface::OpNeumannForceAnalytical::doWork(
108 int side, EntityType type, EntitiesFieldData::EntData &data) {
109 MoFEMFunctionBegin;
110
111 if (data.getIndices().size() == 0)
112 MoFEMFunctionReturnHot(0);
113 const EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
114 if (tRis.find(ent) == tRis.end())
115 MoFEMFunctionReturnHot(0);
116
117 const int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
118 const int nb_row_dofs = data.getIndices().size() / rank;
119
120 nF.resize(data.getIndices().size(), false);
121 nF.clear();
122
123 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
124
125 VectorDouble3 coords(3);
126 VectorDouble3 normal(3);
127 VectorDouble3 force(3);
128
129 for (unsigned int gg = 0; gg != data.getN().size1(); ++gg) {
130
131 // get integration weight and Jacobian of integration point (area of face)
132 double val = getGaussPts()(2, gg);
133 val *= cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
134 if (fe_type == MBTRI)
135 val /= 2;
136 for (int dd = 0; dd != 3; ++dd) {
137 coords[dd] = getCoordsAtGaussPts()(gg, dd);
138 normal[dd] = getNormalsAtGaussPts()(gg, dd);
139 }
140
141 if (analyticalForceOp) {
142 CHKERR analyticalForceOp->getForce(ent, coords, normal, force);
143 for (int rr = 0; rr != 3; ++rr) {
144 cblas_daxpy(nb_row_dofs, val * force[rr], &data.getN()(gg, 0), 1,
145 &nF[rr], rank);
146 }
147 } else {
148 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "No force to apply");
149 }
150 }
151
152 // Scale force using user defined scaling operator
153 CHKERR MethodForForceScaling::applyScale(getFEMethod(), methodsOp, nF);
154
155 auto get_f = [&]() {
156 if (F == PETSC_NULL)
157 return getKSPf();
158 return F;
159 };
160
161 // Assemble force into vector
162 CHKERR VecSetValues(get_f(), data, &*nF.data().begin(), ADD_VALUES);
163
164 MoFEMFunctionReturn(0);
165}
166
167NeumannForcesSurface::OpNeumannPressure::OpNeumannPressure(
168 const std::string field_name, Vec _F, bCPressure &data,
169 boost::ptr_vector<MethodForForceScaling> &methods_op, bool ho_geometry)
170 : FaceElementForcesAndSourcesCore::UserDataOperator(
171 field_name, UserDataOperator::OPROW),
172 F(_F), dAta(data), methodsOp(methods_op), hoGeometry(ho_geometry) {}
173
174MoFEMErrorCode NeumannForcesSurface::OpNeumannPressure::doWork(
175 int side, EntityType type, EntitiesFieldData::EntData &data) {
176
177 MoFEMFunctionBegin;
178
179 if (data.getIndices().size() == 0)
180 MoFEMFunctionReturnHot(0);
181
182 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
183 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
184 if (dAta.tRis.find(fe_ent) == dAta.tRis.end())
185 MoFEMFunctionReturnHot(0);
186
187 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
188 int nb_row_dofs = data.getIndices().size() / rank;
189
190 Nf.resize(data.getIndices().size(), false);
191 Nf.clear();
192
193 for (unsigned int gg = 0; gg != data.getN().size1(); ++gg) {
194
195 double val = getGaussPts()(2, gg);
196 if (fe_type == MBTRI)
197 val /= 2;
198
199 for (int rr = 0; rr != rank; ++rr) {
200
201 double force;
202 if (hoGeometry || fe_type == MBQUAD)
203 force = dAta.data.data.value1 * getNormalsAtGaussPts()(gg, rr);
204 else
205 force = dAta.data.data.value1 * getNormal()[rr];
206
207 cblas_daxpy(nb_row_dofs, val * force, &data.getN()(gg, 0), 1, &Nf[rr],
208 rank);
209 }
210 }
211
212 CHKERR MethodForForceScaling::applyScale(getFEMethod(), methodsOp, Nf);
213
214 auto get_f = [&]() {
215 if (F == PETSC_NULL)
216 return getKSPf();
217 return F;
218 };
219
220 CHKERR VecSetValues(get_f(), data, &*Nf.data().begin(), ADD_VALUES);
221
222 MoFEMFunctionReturn(0);
223}
224
225MoFEMErrorCode
226NeumannForcesSurface::OpGetTangent::doWork(int side, EntityType type,
227 EntitiesFieldData::EntData &data) {
228 MoFEMFunctionBegin;
229
230 if (data.getFieldData().size() == 0)
231 PetscFunctionReturn(0);
232
233 ngp = data.getN().size1();
234
235 unsigned int nb_dofs = data.getFieldData().size() / 3;
236 FTensor::Index<'i', 3> i;
237 if (type == MBVERTEX) {
238 dataAtIntegrationPts->tangent1.resize(3, ngp, false);
239 dataAtIntegrationPts->tangent1.clear();
240
241 dataAtIntegrationPts->tangent2.resize(3, ngp, false);
242 dataAtIntegrationPts->tangent2.clear();
243 }
244
245 auto t_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
246 auto t_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
247
248 for (unsigned int gg = 0; gg != ngp; ++gg) {
249 auto t_N = data.getFTensor1DiffN<2>(gg, 0);
250 auto t_dof = data.getFTensor1FieldData<3>();
251
252 for (unsigned int dd = 0; dd != nb_dofs; ++dd) {
253 t_1(i) += t_dof(i) * t_N(0);
254 t_2(i) += t_dof(i) * t_N(1);
255 ++t_dof;
256 ++t_N;
257 }
258 ++t_1;
259 ++t_2;
260 }
261
262 MoFEMFunctionReturn(0);
263}
264
265MoFEMErrorCode NeumannForcesSurface::OpNeumannPressureLhs_dx_dX::doWork(
266 int row_side, int col_side, EntityType row_type, EntityType col_type,
267 EntitiesFieldData::EntData &row_data,
268 EntitiesFieldData::EntData &col_data) {
269 MoFEMFunctionBegin;
270
271 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
272 dAta.tRis.end()) {
273 MoFEMFunctionReturnHot(0);
274 }
275
276 const int row_nb_dofs = row_data.getIndices().size();
277 if (!row_nb_dofs)
278 MoFEMFunctionReturnHot(0);
279 const int col_nb_dofs = col_data.getIndices().size();
280 if (!col_nb_dofs)
281 MoFEMFunctionReturnHot(0);
282 const int nb_gauss_pts = row_data.getN().size1();
283
284 int nb_base_fun_row = row_data.getFieldData().size() / 3;
285 int nb_base_fun_col = col_data.getFieldData().size() / 3;
286
287 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
288 NN.clear();
289
290 FTensor::Index<'i', 3> i;
291 FTensor::Index<'j', 3> j;
292 FTensor::Index<'k', 3> k;
293
294 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
295 return FTensor::Tensor2<double *, 3, 3>(
296 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
297 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
298 &m(r + 2, c + 2));
299 };
300
301 auto make_vec_der = [](auto t_N, auto t_1, auto t_2) {
302 FTensor::Index<'i', 3> i;
303 FTensor::Index<'j', 3> j;
304 FTensor::Index<'k', 3> k;
305 FTensor::Tensor2<double, 3, 3> t_n;
306 t_n(i, j) = 0;
307 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
308 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
309 return t_n;
310 };
311
312 double lambda = 1;
313 if (auto arc_dof = dataAtIntegrationPts->arcLengthDof.lock()) {
314 lambda = arc_dof->getFieldData();
315 }
316
317 auto t_w = getFTensor0IntegrationWeight();
318 auto t_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
319 auto t_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
320 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
321
322 double val = 0.5 * t_w * lambda;
323
324 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
325
326 int bbc = 0;
327 for (; bbc != nb_base_fun_col; ++bbc) {
328 auto t_base = row_data.getFTensor0N(gg, 0);
329
330 int bbr = 0;
331 for (; bbr != nb_base_fun_row; ++bbr) {
332
333 auto t_d_n = make_vec_der(t_N, t_1, t_2);
334
335 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
336
337 // TODO: handle hoGeometry
338 t_assemble(i, k) += val * dAta.data.data.value1 * t_base * t_d_n(i, k);
339
340 ++t_base;
341 }
342 ++t_N;
343 }
344 ++t_w;
345 ++t_1;
346 ++t_2;
347 }
348
349 // get pointer to first global index on row
350 const int *row_indices = &*row_data.getIndices().data().begin();
351 // get pointer to first global index on column
352 const int *col_indices = &*col_data.getIndices().data().begin();
353
354 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
355 : getFEMethod()->snes_B;
356
357 CHKERR MatSetValues(B, row_nb_dofs, row_indices, col_nb_dofs, col_indices,
358 &*NN.data().begin(), ADD_VALUES);
359
360 MoFEMFunctionReturn(0);
361}
362
363MoFEMErrorCode NeumannForcesSurface::OpNeumannSurfaceForceLhs_dx_dX::doWork(
364 int row_side, int col_side, EntityType row_type, EntityType col_type,
365 EntitiesFieldData::EntData &row_data,
366 EntitiesFieldData::EntData &col_data) {
367 MoFEMFunctionBegin;
368
369 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
370 dAta.tRis.end()) {
371 MoFEMFunctionReturnHot(0);
372 }
373
374 const int row_nb_dofs = row_data.getIndices().size();
375 if (!row_nb_dofs)
376 MoFEMFunctionReturnHot(0);
377 const int col_nb_dofs = col_data.getIndices().size();
378 if (!col_nb_dofs)
379 MoFEMFunctionReturnHot(0);
380 const int nb_gauss_pts = row_data.getN().size1();
381
382 int nb_base_fun_row = row_data.getFieldData().size() / 3;
383 int nb_base_fun_col = col_data.getFieldData().size() / 3;
384
385 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
386 NN.clear();
387
388 FTensor::Index<'i', 3> i;
389 FTensor::Index<'j', 3> j;
390 FTensor::Index<'k', 3> k;
391
392 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
393 return FTensor::Tensor2<double *, 3, 3>(
394 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
395 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
396 &m(r + 2, c + 2));
397 };
398
399 auto make_vec_der = [](auto t_N, auto t_1, auto t_2) {
400 FTensor::Index<'i', 3> i;
401 FTensor::Index<'j', 3> j;
402 FTensor::Index<'k', 3> k;
403 FTensor::Tensor2<double, 3, 3> t_n;
404 t_n(i, j) = 0;
405 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
406 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
407 return t_n;
408 };
409
410 double lambda = 1;
411 if (auto arc_dof = dataAtIntegrationPts->arcLengthDof.lock()) {
412 lambda = arc_dof->getFieldData();
413 }
414
415 auto t_w = getFTensor0IntegrationWeight();
416 auto t_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
417 auto t_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
418 auto t_normal = getFTensor1NormalsAtGaussPts();
419 FTensor1 t_force(dAta.data.data.value3, dAta.data.data.value4,
420 dAta.data.data.value5);
421 t_force(i) *= dAta.data.data.value1;
422
423 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
424
425 double nrm2 = std::sqrt(t_normal(i) * t_normal(i));
426
427 double val = 0.5 * t_w * lambda / nrm2;
428
429 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
430
431 int bbc = 0;
432 for (; bbc != nb_base_fun_col; ++bbc) {
433 auto t_base = row_data.getFTensor0N(gg, 0);
434
435 int bbr = 0;
436 for (; bbr != nb_base_fun_row; ++bbr) {
437
438 auto t_d_n = make_vec_der(t_N, t_1, t_2);
439
440 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
441
442 // TODO: handle hoGeometry
443 t_assemble(i, k) +=
444 val * t_force(i) * t_base * t_d_n(j, k) * t_normal(j);
445
446 ++t_base;
447 }
448 ++t_N;
449 }
450 ++t_w;
451 ++t_1;
452 ++t_2;
453 ++t_normal;
454 }
455
456 // get pointer to first global index on row
457 const int *row_indices = &*row_data.getIndices().data().begin();
458 // get pointer to first global index on column
459 const int *col_indices = &*col_data.getIndices().data().begin();
460
461 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
462 : getFEMethod()->snes_B;
463
464 CHKERR MatSetValues(B, row_nb_dofs, row_indices, col_nb_dofs, col_indices,
465 &*NN.data().begin(), ADD_VALUES);
466
467 MoFEMFunctionReturn(0);
468}
469
470MoFEMErrorCode NeumannForcesSurface::OpCalculateDeformation::doWork(
471 int side, EntityType type, EntitiesFieldData::EntData &row_data) {
472
473 MoFEMFunctionBegin;
474 FTensor::Index<'i', 3> i;
475 FTensor::Index<'j', 3> j;
476 FTensor::Index<'k', 3> k;
477
478 // get number of integration points
479 const int nb_integration_pts = getGaussPts().size2();
480 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hMat);
481 auto t_H = getFTensor2FromMat<3, 3>(dataAtPts->HMat);
482
483 dataAtPts->detHVec.resize(nb_integration_pts, false);
484 dataAtPts->invHMat.resize(9, nb_integration_pts, false);
485 dataAtPts->FMat.resize(9, nb_integration_pts, false);
486
487 auto t_detH = getFTensor0FromVec(dataAtPts->detHVec);
488 auto t_invH = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
489 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
490
491 for (int gg = 0; gg != nb_integration_pts; ++gg) {
492 CHKERR determinantTensor3by3(t_H, t_detH);
493 CHKERR invertTensor3by3(t_H, t_detH, t_invH);
494 t_F(i, j) = t_h(i, k) * t_invH(k, j);
495
496 ++t_h;
497 ++t_H;
498 ++t_detH;
499 ++t_invH;
500 ++t_F;
501 }
502
503 MoFEMFunctionReturn(0);
504}
505
506MoFEMErrorCode NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX::doWork(
507 int side, EntityType type, EntitiesFieldData::EntData &row_data) {
508
509 MoFEMFunctionBegin;
510
511 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
512 dAta.tRis.end()) {
513 MoFEMFunctionReturnHot(0);
514 }
515
516 // get number of dofs on row
517 nbRows = row_data.getIndices().size();
518 // if no dofs on row, exit that work, nothing to do here
519 if (!nbRows)
520 MoFEMFunctionReturnHot(0);
521
522 nF.resize(nbRows, false);
523 nF.clear();
524
525 // get number of integration points
526 nbIntegrationPts = getGaussPts().size2();
527
528 // integrate local matrix for entity block
529 CHKERR iNtegrate(row_data);
530
531 // assemble local matrix
532 CHKERR aSsemble(row_data);
533
534 MoFEMFunctionReturn(0);
535}
536
537MoFEMErrorCode NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX::iNtegrate(
538 EntData &data) {
539 MoFEMFunctionBegin;
540
541 auto get_tensor1 = [](VectorDouble &v, const int r) {
542 return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
543 &v(r + 0), &v(r + 1), &v(r + 2));
544 };
545
546 FTensor::Index<'i', 3> i;
547 FTensor::Index<'j', 3> j;
548
549 CHKERR loopSideVolumes(sideFeName, *sideFe);
550
551 // get integration weights
552 auto t_w = getFTensor0IntegrationWeight();
553
554 auto t_normal = getFTensor1NormalsAtGaussPts();
555
556 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
557
558 // iterate over integration points
559 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
560
561 FTensor::Tensor0<double *> t_base(&data.getN()(gg, 0));
562
563 double a = -0.5 * t_w * dAta.data.data.value1;
564 auto t_nf = get_tensor1(nF, 0);
565
566 int rr = 0;
567 for (; rr != nbRows / 3; ++rr) {
568
569 t_nf(i) += a * t_base * t_F(j, i) * t_normal(j);
570
571 ++t_base;
572 ++t_nf;
573 }
574
575 ++t_w;
576 ++t_F;
577 ++t_normal;
578 }
579
580 MoFEMFunctionReturn(0);
581}
582
583MoFEMErrorCode NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX::aSsemble(
584 EntData &row_data) {
585 MoFEMFunctionBegin;
586
587 // get pointer to first global index on row
588 const int *row_indices = &*row_data.getIndices().data().begin();
589
590 auto &data = *dataAtPts;
591
592 rowIndices.resize(nbRows, false);
593 noalias(rowIndices) = row_data.getIndices();
594 row_indices = &rowIndices[0];
595 VectorDofs &dofs = row_data.getFieldDofs();
596 VectorDofs::iterator dit = dofs.begin();
597 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
598 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
599 data.forcesOnlyOnEntitiesRow.end()) {
600 rowIndices[ii] = -1;
601 }
602 }
603
604 auto get_f = [&]() {
605 if (F == PETSC_NULL)
606 return getKSPf();
607 return F;
608 };
609
610 auto vec_assemble = [&](Vec my_f) {
611 MoFEMFunctionBegin;
612 CHKERR VecSetOption(my_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
613 CHKERR VecSetValues(my_f, nbRows, row_indices, &*nF.data().begin(),
614 ADD_VALUES);
615 MoFEMFunctionReturn(0);
616 };
617
618 // assemble local matrix
619 CHKERR vec_assemble(get_f());
620
621 MoFEMFunctionReturn(0);
622}
623
624MoFEMErrorCode
625NeumannForcesSurface::OpNeumannSurfaceForceMaterialRhs_dX::doWork(
626 int side, EntityType type, EntitiesFieldData::EntData &row_data) {
627
628 MoFEMFunctionBegin;
629
630 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
631 dAta.tRis.end()) {
632 MoFEMFunctionReturnHot(0);
633 }
634
635 // get number of dofs on row
636 nbRows = row_data.getIndices().size();
637 // if no dofs on row, exit that work, nothing to do here
638 if (!nbRows)
639 MoFEMFunctionReturnHot(0);
640
641 nF.resize(nbRows, false);
642 nF.clear();
643
644 // get number of integration points
645 nbIntegrationPts = getGaussPts().size2();
646
647 // integrate local matrix for entity block
648 CHKERR iNtegrate(row_data);
649
650 // assemble local matrix
651 CHKERR aSsemble(row_data);
652
653 MoFEMFunctionReturn(0);
654}
655
656MoFEMErrorCode
657NeumannForcesSurface::OpNeumannSurfaceForceMaterialRhs_dX::iNtegrate(
658 EntData &data) {
659 MoFEMFunctionBegin;
660
661 auto get_tensor1 = [](VectorDouble &v, const int r) {
662 return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
663 &v(r + 0), &v(r + 1), &v(r + 2));
664 };
665
666 FTensor::Index<'i', 3> i;
667 FTensor::Index<'j', 3> j;
668
669 CHKERR loopSideVolumes(sideFeName, *sideFe);
670
671 // get integration weights
672 auto t_w = getFTensor0IntegrationWeight();
673
674 auto t_normal = getFTensor1NormalsAtGaussPts();
675
676 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
677
678 FTensor1 t_force(dAta.data.data.value3, dAta.data.data.value4,
679 dAta.data.data.value5);
680 t_force(i) *= dAta.data.data.value1;
681 // iterate over integration points
682 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
683
684 FTensor::Tensor0<double *> t_base(&data.getN()(gg, 0));
685
686 double nrm2 = std::sqrt(t_normal(i) * t_normal(i));
687
688 double a = -0.5 * t_w * nrm2;
689 auto t_nf = get_tensor1(nF, 0);
690
691 int rr = 0;
692 for (; rr != nbRows / 3; ++rr) {
693
694 t_nf(i) += a * t_base * t_F(j, i) * t_force(j);
695
696 ++t_base;
697 ++t_nf;
698 }
699
700 ++t_w;
701 ++t_F;
702 ++t_normal;
703 }
704
705 MoFEMFunctionReturn(0);
706}
707
708MoFEMErrorCode
709NeumannForcesSurface::OpNeumannSurfaceForceMaterialRhs_dX::aSsemble(
710 EntData &row_data) {
711 MoFEMFunctionBegin;
712
713 // get pointer to first global index on row
714 const int *row_indices = &*row_data.getIndices().data().begin();
715
716 auto &data = *dataAtPts;
717
718 rowIndices.resize(nbRows, false);
719 noalias(rowIndices) = row_data.getIndices();
720 row_indices = &rowIndices[0];
721 VectorDofs &dofs = row_data.getFieldDofs();
722 VectorDofs::iterator dit = dofs.begin();
723 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
724 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
725 data.forcesOnlyOnEntitiesRow.end()) {
726 rowIndices[ii] = -1;
727 }
728 }
729
730 auto get_f = [&]() {
731 if (F == PETSC_NULL)
732 return getKSPf();
733 return F;
734 };
735
736 auto vec_assemble = [&](Vec my_f) {
737 MoFEMFunctionBegin;
738 CHKERR VecSetOption(my_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
739 CHKERR VecSetValues(my_f, nbRows, row_indices, &*nF.data().begin(),
740 ADD_VALUES);
741 MoFEMFunctionReturn(0);
742 };
743
744 // assemble local matrix
745 CHKERR vec_assemble(get_f());
746
747 MoFEMFunctionReturn(0);
748}
749
750MoFEMErrorCode NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dX::doWork(
751 int row_side, int col_side, EntityType row_type, EntityType col_type,
752 EntitiesFieldData::EntData &row_data,
753 EntitiesFieldData::EntData &col_data) {
754
755 MoFEMFunctionBegin;
756
757 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
758 dAta.tRis.end()) {
759 MoFEMFunctionReturnHot(0);
760 }
761
762 row_nb_dofs = row_data.getIndices().size();
763 if (!row_nb_dofs)
764 MoFEMFunctionReturnHot(0);
765 col_nb_dofs = col_data.getIndices().size();
766 if (!col_nb_dofs)
767 MoFEMFunctionReturnHot(0);
768 nb_gauss_pts = row_data.getN().size1();
769
770 nb_base_fun_row = row_data.getFieldData().size() / 3;
771 nb_base_fun_col = col_data.getFieldData().size() / 3;
772
773 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
774 NN.clear();
775
776 diagonal_block = (row_type == col_type) && (row_side == col_side);
777
778 if (col_type == MBVERTEX) {
779 dataAtPts->faceRowData = &row_data;
780 CHKERR loopSideVolumes(sideFeName, *sideFe);
781 }
782
783 // integrate local matrix for entity block
784 CHKERR iNtegrate(row_data, col_data);
785
786 // assemble local matrix
787 CHKERR aSsemble(row_data, col_data);
788
789 MoFEMFunctionReturn(0);
790}
791
792MoFEMErrorCode
793NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dX::iNtegrate(
794 EntData &row_data, EntData &col_data) {
795
796 MoFEMFunctionBegin;
797
798 FTensor::Index<'i', 3> i;
799 FTensor::Index<'j', 3> j;
800 FTensor::Index<'k', 3> k;
801
802 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
803 return FTensor::Tensor2<double *, 3, 3>(
804 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
805 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
806 &m(r + 2, c + 2));
807 };
808
809 auto make_vec_der = [](auto t_N, auto t_1, auto t_2) {
810 FTensor::Index<'i', 3> i;
811 FTensor::Index<'j', 3> j;
812 FTensor::Index<'k', 3> k;
813 FTensor::Tensor2<double, 3, 3> t_n;
814 t_n(i, j) = 0;
815 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
816 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
817 return t_n;
818 };
819
820 dataAtPts->faceRowData = nullptr;
821 CHKERR loopSideVolumes(sideFeName, *sideFe);
822
823 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
824
825 double lambda = 1;
826 if (auto arc_dof = dataAtPts->arcLengthDof.lock()) {
827 lambda = arc_dof->getFieldData();
828 }
829
830 auto t_w = getFTensor0IntegrationWeight();
831 auto t_1 = getFTensor1FromMat<3>(dataAtPts->tangent1);
832 auto t_2 = getFTensor1FromMat<3>(dataAtPts->tangent2);
833
834 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
835
836 double val = 0.5 * t_w * lambda;
837
838 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
839
840 int bbc = 0;
841 for (; bbc != nb_base_fun_col; ++bbc) {
842
843 auto t_base = row_data.getFTensor0N(gg, 0);
844
845 int bbr = 0;
846 for (; bbr != nb_base_fun_row; ++bbr) {
847
848 auto t_d_n = make_vec_der(t_N, t_1, t_2);
849
850 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
851 // TODO: handle hoGeometry
852
853 t_assemble(i, k) -=
854 val * dAta.data.data.value1 * t_base * t_F(j, i) * t_d_n(j, k);
855
856 ++t_base;
857 }
858 ++t_N;
859 }
860 ++t_F;
861 ++t_w;
862 ++t_1;
863 ++t_2;
864 }
865
866 MoFEMFunctionReturn(0);
867}
868
869MoFEMErrorCode NeumannForcesSurface::OpNeumannSurfaceForceMaterialLhs::aSsemble(
870 EntData &row_data, EntData &col_data) {
871
872 MoFEMFunctionBegin;
873
874 // get pointer to first global index on row
875 const int *row_indices = &*row_data.getIndices().data().begin();
876 // get pointer to first global index on column
877 const int *col_indices = &*col_data.getIndices().data().begin();
878
879 auto &data = *dataAtPts;
880
881 rowIndices.resize(row_nb_dofs, false);
882 noalias(rowIndices) = row_data.getIndices();
883 row_indices = &rowIndices[0];
884 VectorDofs &dofs = row_data.getFieldDofs();
885 VectorDofs::iterator dit = dofs.begin();
886 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
887 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
888 data.forcesOnlyOnEntitiesRow.end()) {
889 rowIndices[ii] = -1;
890 }
891 }
892
893 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
894 : getFEMethod()->snes_B;
895 // assemble local matrix
896 CHKERR MatSetValues(B, row_nb_dofs, row_indices, col_nb_dofs, col_indices,
897 &*NN.data().begin(), ADD_VALUES);
898
899 MoFEMFunctionReturn(0);
900}
901
902MoFEMErrorCode
903NeumannForcesSurface::OpNeumannSurfaceForceMaterialLhs_dX_dX::doWork(
904 int row_side, int col_side, EntityType row_type, EntityType col_type,
905 EntitiesFieldData::EntData &row_data,
906 EntitiesFieldData::EntData &col_data) {
907
908 MoFEMFunctionBegin;
909
910 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
911 dAta.tRis.end()) {
912 MoFEMFunctionReturnHot(0);
913 }
914
915 row_nb_dofs = row_data.getIndices().size();
916 if (!row_nb_dofs)
917 MoFEMFunctionReturnHot(0);
918 col_nb_dofs = col_data.getIndices().size();
919 if (!col_nb_dofs)
920 MoFEMFunctionReturnHot(0);
921 nb_gauss_pts = row_data.getN().size1();
922
923 nb_base_fun_row = row_data.getFieldData().size() / 3;
924 nb_base_fun_col = col_data.getFieldData().size() / 3;
925
926 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
927 NN.clear();
928
929 diagonal_block = (row_type == col_type) && (row_side == col_side);
930
931 if (col_type == MBVERTEX) {
932 dataAtPts->faceRowData = &row_data;
933 CHKERR loopSideVolumes(sideFeName, *sideFe);
934 }
935
936 // integrate local matrix for entity block
937 CHKERR iNtegrate(row_data, col_data);
938
939 // assemble local matrix
940 CHKERR aSsemble(row_data, col_data);
941
942 MoFEMFunctionReturn(0);
943}
944
945MoFEMErrorCode
946NeumannForcesSurface::OpNeumannSurfaceForceMaterialLhs_dX_dX::iNtegrate(
947 EntData &row_data, EntData &col_data) {
948
949 MoFEMFunctionBegin;
950
951 FTensor::Index<'i', 3> i;
952 FTensor::Index<'j', 3> j;
953 FTensor::Index<'k', 3> k;
954 FTensor::Index<'l', 3> l;
955
956 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
957 return FTensor::Tensor2<double *, 3, 3>(
958 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
959 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
960 &m(r + 2, c + 2));
961 };
962
963 auto make_vec_der = [](auto t_N, auto t_1, auto t_2) {
964 FTensor::Index<'i', 3> i;
965 FTensor::Index<'j', 3> j;
966 FTensor::Index<'k', 3> k;
967 FTensor::Tensor2<double, 3, 3> t_n;
968 t_n(i, j) = 0;
969 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
970 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
971 return t_n;
972 };
973
974 dataAtPts->faceRowData = nullptr;
975 CHKERR loopSideVolumes(sideFeName, *sideFe);
976
977 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
978
979 double lambda = 1;
980 if (auto arc_dof = dataAtPts->arcLengthDof.lock()) {
981 lambda = arc_dof->getFieldData();
982 }
983
984 auto t_w = getFTensor0IntegrationWeight();
985 auto t_1 = getFTensor1FromMat<3>(dataAtPts->tangent1);
986 auto t_2 = getFTensor1FromMat<3>(dataAtPts->tangent2);
987
988 auto t_normal = getFTensor1NormalsAtGaussPts();
989 FTensor1 t_force(dAta.data.data.value3, dAta.data.data.value4,
990 dAta.data.data.value5);
991 t_force(i) *= dAta.data.data.value1;
992 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
993
994 double nrm2 = std::sqrt(t_normal(i) * t_normal(i));
995
996 double val = 0.5 * t_w * lambda / nrm2;
997
998 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
999
1000 int bbc = 0;
1001 for (; bbc != nb_base_fun_col; ++bbc) {
1002
1003 auto t_base = row_data.getFTensor0N(gg, 0);
1004
1005 int bbr = 0;
1006 for (; bbr != nb_base_fun_row; ++bbr) {
1007
1008 auto t_d_n = make_vec_der(t_N, t_1, t_2);
1009
1010 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
1011 // TODO: handle hoGeometry
1012
1013 t_assemble(i, k) -=
1014 val * t_base * t_force(j) * t_F(j, i) * (t_d_n(l, k) * t_normal(l));
1015
1016 ++t_base;
1017 }
1018 ++t_N;
1019 }
1020 ++t_F;
1021 ++t_w;
1022 ++t_1;
1023 ++t_2;
1024 ++t_normal;
1025 }
1026
1027 MoFEMFunctionReturn(0);
1028}
1029
1030MoFEMErrorCode NeumannForcesSurface::OpNeumannPressureMaterialLhs::aSsemble(
1031 EntData &row_data, EntData &col_data) {
1032
1033 MoFEMFunctionBegin;
1034
1035 // get pointer to first global index on row
1036 const int *row_indices = &*row_data.getIndices().data().begin();
1037 // get pointer to first global index on column
1038 const int *col_indices = &*col_data.getIndices().data().begin();
1039
1040 auto &data = *dataAtPts;
1041
1042 rowIndices.resize(row_nb_dofs, false);
1043 noalias(rowIndices) = row_data.getIndices();
1044 row_indices = &rowIndices[0];
1045 VectorDofs &dofs = row_data.getFieldDofs();
1046 VectorDofs::iterator dit = dofs.begin();
1047 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
1048 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
1049 data.forcesOnlyOnEntitiesRow.end()) {
1050 rowIndices[ii] = -1;
1051 }
1052 }
1053
1054 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
1055 : getFEMethod()->snes_B;
1056 // assemble local matrix
1057 CHKERR MatSetValues(B, row_nb_dofs, row_indices, col_nb_dofs, col_indices,
1058 &*NN.data().begin(), ADD_VALUES);
1059
1060 MoFEMFunctionReturn(0);
1061}
1062
1063MoFEMErrorCode NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dx::doWork(
1064 int row_side, int col_side, EntityType row_type, EntityType col_type,
1065 EntitiesFieldData::EntData &row_data,
1066 EntitiesFieldData::EntData &col_data) {
1067
1068 MoFEMFunctionBegin;
1069
1070 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1071 dAta.tRis.end()) {
1072 MoFEMFunctionReturnHot(0);
1073 }
1074
1075 if (col_type != MBVERTEX)
1076 MoFEMFunctionReturnHot(0);
1077
1078 dataAtPts->faceRowData = &row_data;
1079 CHKERR loopSideVolumes(sideFeName, *sideFe);
1080
1081 MoFEMFunctionReturn(0);
1082}
1083
1084MoFEMErrorCode
1085NeumannForcesSurface::OpNeumannSurfaceForceMaterialLhs_dX_dx::doWork(
1086 int row_side, int col_side, EntityType row_type, EntityType col_type,
1087 EntitiesFieldData::EntData &row_data,
1088 EntitiesFieldData::EntData &col_data) {
1089
1090 MoFEMFunctionBegin;
1091
1092 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1093 dAta.tRis.end()) {
1094 MoFEMFunctionReturnHot(0);
1095 }
1096
1097 if (col_type != MBVERTEX)
1098 MoFEMFunctionReturnHot(0);
1099
1100 dataAtPts->faceRowData = &row_data;
1101 CHKERR loopSideVolumes(sideFeName, *sideFe);
1102
1103 MoFEMFunctionReturn(0);
1104}
1105
1106MoFEMErrorCode
1107NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs::doWork(
1108 int row_side, int col_side, EntityType row_type, EntityType col_type,
1109 EntitiesFieldData::EntData &row_data,
1110 EntitiesFieldData::EntData &col_data) {
1111
1112 MoFEMFunctionBegin;
1113
1114 if (dataAtPts->faceRowData == nullptr)
1115 MoFEMFunctionReturnHot(0);
1116
1117 if (row_type != MBVERTEX)
1118 MoFEMFunctionReturnHot(0);
1119
1120 row_nb_dofs = dataAtPts->faceRowData->getIndices().size();
1121 if (!row_nb_dofs)
1122 MoFEMFunctionReturnHot(0);
1123 col_nb_dofs = col_data.getIndices().size();
1124 if (!col_nb_dofs)
1125 MoFEMFunctionReturnHot(0);
1126
1127 nb_gauss_pts = dataAtPts->faceRowData->getN().size1();
1128
1129 nb_base_fun_row = dataAtPts->faceRowData->getFieldData().size() / 3;
1130 nb_base_fun_col = col_data.getFieldData().size() / 3;
1131
1132 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
1133 NN.clear();
1134
1135 // integrate local matrix for entity block
1136 CHKERR iNtegrate(*(dataAtPts->faceRowData), col_data);
1137
1138 // assemble local matrix
1139 CHKERR aSsemble(*(dataAtPts->faceRowData), col_data);
1140
1141 MoFEMFunctionReturn(0);
1142}
1143
1144MoFEMErrorCode
1145NeumannForcesSurface::OpNeumannSurfaceForceMaterialVolOnSideLhs::doWork(
1146 int row_side, int col_side, EntityType row_type, EntityType col_type,
1147 EntitiesFieldData::EntData &row_data,
1148 EntitiesFieldData::EntData &col_data) {
1149
1150 MoFEMFunctionBegin;
1151
1152 if (dataAtPts->faceRowData == nullptr)
1153 MoFEMFunctionReturnHot(0);
1154
1155 if (row_type != MBVERTEX)
1156 MoFEMFunctionReturnHot(0);
1157
1158 row_nb_dofs = dataAtPts->faceRowData->getIndices().size();
1159 if (!row_nb_dofs)
1160 MoFEMFunctionReturnHot(0);
1161 col_nb_dofs = col_data.getIndices().size();
1162 if (!col_nb_dofs)
1163 MoFEMFunctionReturnHot(0);
1164
1165 nb_gauss_pts = dataAtPts->faceRowData->getN().size1();
1166
1167 nb_base_fun_row = dataAtPts->faceRowData->getFieldData().size() / 3;
1168 nb_base_fun_col = col_data.getFieldData().size() / 3;
1169
1170 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
1171 NN.clear();
1172
1173 // integrate local matrix for entity block
1174 CHKERR iNtegrate(*(dataAtPts->faceRowData), col_data);
1175
1176 // assemble local matrix
1177 CHKERR aSsemble(*(dataAtPts->faceRowData), col_data);
1178
1179 MoFEMFunctionReturn(0);
1180}
1181
1182MoFEMErrorCode
1183NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs_dX_dx::iNtegrate(
1184 EntData &row_data, EntData &col_data) {
1185
1186 MoFEMFunctionBegin;
1187
1188 FTensor::Index<'i', 3> i;
1189 FTensor::Index<'j', 3> j;
1190 FTensor::Index<'k', 3> k;
1191
1192 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1193 return FTensor::Tensor2<double *, 3, 3>(
1194 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1195 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1196 &m(r + 2, c + 2));
1197 };
1198
1199 auto t_w = getFTensor0IntegrationWeight();
1200
1201 auto t_normal = getFTensor1NormalsAtGaussPts();
1202
1203 auto t_inv_H = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
1204
1205 double lambda = 1;
1206 if (auto arc_dof = dataAtPts->arcLengthDof.lock()) {
1207 lambda = arc_dof->getFieldData();
1208 }
1209
1210 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1211
1212 double a = -0.5 * t_w * dAta.data.data.value1 * lambda;
1213
1214 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1215
1216 int bbc = 0;
1217 for (; bbc != nb_base_fun_col; ++bbc) {
1218
1219 FTensor::Tensor0<double *> t_row_base(&row_data.getN()(gg, 0));
1220
1221 int bbr = 0;
1222 for (; bbr != nb_base_fun_row; ++bbr) {
1223
1224 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
1225 // TODO: handle hoGeometry
1226
1227 t_assemble(i, j) +=
1228 a * t_row_base * t_inv_H(k, i) * t_col_diff_base(k) * t_normal(j);
1229
1230 ++t_row_base;
1231 }
1232 ++t_col_diff_base;
1233 }
1234 ++t_w;
1235 ++t_normal;
1236 ++t_inv_H;
1237 }
1238
1239 MoFEMFunctionReturn(0);
1240}
1241
1242MoFEMErrorCode
1243NeumannForcesSurface::OpNeumannSurfaceForceMaterialVolOnSideLhs_dX_dx::
1244 iNtegrate(EntData &row_data, EntData &col_data) {
1245
1246 MoFEMFunctionBegin;
1247
1248 FTensor::Index<'i', 3> i;
1249 FTensor::Index<'j', 3> j;
1250 FTensor::Index<'k', 3> k;
1251
1252 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1253 return FTensor::Tensor2<double *, 3, 3>(
1254 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1255 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1256 &m(r + 2, c + 2));
1257 };
1258
1259 auto t_w = getFTensor0IntegrationWeight();
1260
1261 auto t_normal = getFTensor1NormalsAtGaussPts();
1262
1263 auto t_inv_H = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
1264
1265 double lambda = 1;
1266 if (auto arc_dof = dataAtPts->arcLengthDof.lock()) {
1267 lambda = arc_dof->getFieldData();
1268 }
1269
1270 FTensor1 t_force(dAta.data.data.value3, dAta.data.data.value4,
1271 dAta.data.data.value5);
1272 t_force(i) *= dAta.data.data.value1;
1273 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1274
1275 double nrm2 = std::sqrt(t_normal(i) * t_normal(i));
1276
1277 double a = -0.5 * t_w * nrm2 * lambda;
1278
1279 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1280
1281 int bbc = 0;
1282 for (; bbc != nb_base_fun_col; ++bbc) {
1283
1284 FTensor::Tensor0<double *> t_row_base(&row_data.getN()(gg, 0));
1285
1286 int bbr = 0;
1287 for (; bbr != nb_base_fun_row; ++bbr) {
1288
1289 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
1290 // TODO: handle hoGeometry
1291
1292 t_assemble(i, j) +=
1293 a * t_row_base * t_inv_H(k, i) * t_col_diff_base(k) * t_force(j);
1294
1295 ++t_row_base;
1296 }
1297 ++t_col_diff_base;
1298 }
1299 ++t_w;
1300 ++t_normal;
1301 ++t_inv_H;
1302 }
1303
1304 MoFEMFunctionReturn(0);
1305}
1306
1307MoFEMErrorCode
1308NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs::aSsemble(
1309 EntData &row_data, EntData &col_data) {
1310
1311 MoFEMFunctionBegin;
1312
1313 // get pointer to first global index on row
1314 const int *row_indices = &*row_data.getIndices().data().begin();
1315 // get pointer to first global index on column
1316 const int *col_indices = &*col_data.getIndices().data().begin();
1317
1318 auto &data = *dataAtPts;
1319 rowIndices.resize(row_nb_dofs, false);
1320 noalias(rowIndices) = row_data.getIndices();
1321 row_indices = &rowIndices[0];
1322 VectorDofs &dofs = row_data.getFieldDofs();
1323 VectorDofs::iterator dit = dofs.begin();
1324 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
1325 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
1326 data.forcesOnlyOnEntitiesRow.end()) {
1327 rowIndices[ii] = -1;
1328 }
1329 }
1330
1331 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
1332 : getFEMethod()->snes_B;
1333 // assemble local matrix
1334 CHKERR MatSetValues(B, row_nb_dofs, row_indices, col_nb_dofs, col_indices,
1335 &*NN.data().begin(), ADD_VALUES);
1336
1337 MoFEMFunctionReturn(0);
1338}
1339
1340MoFEMErrorCode
1341NeumannForcesSurface::OpNeumannSurfaceForceMaterialVolOnSideLhs::aSsemble(
1342 EntData &row_data, EntData &col_data) {
1343
1344 MoFEMFunctionBegin;
1345
1346 // get pointer to first global index on row
1347 const int *row_indices = &*row_data.getIndices().data().begin();
1348 // get pointer to first global index on column
1349 const int *col_indices = &*col_data.getIndices().data().begin();
1350
1351 auto &data = *dataAtPts;
1352
1353 rowIndices.resize(row_nb_dofs, false);
1354 noalias(rowIndices) = row_data.getIndices();
1355 row_indices = &rowIndices[0];
1356 VectorDofs &dofs = row_data.getFieldDofs();
1357 VectorDofs::iterator dit = dofs.begin();
1358 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
1359 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
1360 data.forcesOnlyOnEntitiesRow.end()) {
1361 rowIndices[ii] = -1;
1362 }
1363 }
1364
1365 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
1366 : getFEMethod()->snes_B;
1367
1368 // assemble local matrix
1369 CHKERR MatSetValues(B, row_nb_dofs, row_indices, col_nb_dofs, col_indices,
1370 &*NN.data().begin(), ADD_VALUES);
1371
1372 MoFEMFunctionReturn(0);
1373}
1374
1375MoFEMErrorCode
1376NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs_dX_dX::iNtegrate(
1377 EntData &row_data, EntData &col_data) {
1378
1379 MoFEMFunctionBegin;
1380
1381 FTensor::Index<'i', 3> i;
1382 FTensor::Index<'j', 3> j;
1383 FTensor::Index<'k', 3> k;
1384 FTensor::Index<'l', 3> l;
1385 FTensor::Index<'m', 3> m;
1386
1387 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1388 return FTensor::Tensor2<double *, 3, 3>(
1389 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1390 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1391 &m(r + 2, c + 2));
1392 };
1393
1394 auto t_w = getFTensor0IntegrationWeight();
1395
1396 auto t_normal = getFTensor1NormalsAtGaussPts();
1397
1398 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hMat);
1399 auto t_inv_H = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
1400
1401 FTensor::Tensor2<double, 3, 3> t_d;
1402
1403 double lambda = 1;
1404 if (auto arc_dof = dataAtPts->arcLengthDof.lock()) {
1405 lambda = arc_dof->getFieldData();
1406 }
1407
1408 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1409
1410 double a = -0.5 * t_w * dAta.data.data.value1 * lambda;
1411
1412 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1413
1414 int bbc = 0;
1415 for (; bbc != nb_base_fun_col; ++bbc) {
1416
1417 FTensor::Tensor0<double *> t_row_base(&row_data.getN()(gg, 0));
1418
1419 int bbr = 0;
1420 for (; bbr != nb_base_fun_row; ++bbr) {
1421
1422 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
1423
1424 // TODO: handle hoGeometry
1425
1426 t_assemble(i, j) -= a * t_row_base * t_inv_H(l, j) *
1427 t_col_diff_base(m) * t_inv_H(m, i) * t_h(k, l) *
1428 t_normal(k);
1429
1430 ++t_row_base;
1431 }
1432 ++t_col_diff_base;
1433 }
1434 ++t_w;
1435 ++t_h;
1436 ++t_inv_H;
1437 ++t_normal;
1438 }
1439
1440 MoFEMFunctionReturn(0);
1441}
1442
1443MoFEMErrorCode
1444NeumannForcesSurface::OpNeumannSurfaceForceMaterialVolOnSideLhs_dX_dX::
1445 iNtegrate(EntData &row_data, EntData &col_data) {
1446
1447 MoFEMFunctionBegin;
1448
1449 FTensor::Index<'i', 3> i;
1450 FTensor::Index<'j', 3> j;
1451 FTensor::Index<'k', 3> k;
1452 FTensor::Index<'l', 3> l;
1453 FTensor::Index<'m', 3> m;
1454
1455 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1456 return FTensor::Tensor2<double *, 3, 3>(
1457 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1458 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1459 &m(r + 2, c + 2));
1460 };
1461
1462 auto t_w = getFTensor0IntegrationWeight();
1463
1464 auto t_normal = getFTensor1NormalsAtGaussPts();
1465
1466 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hMat);
1467 auto t_inv_H = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
1468
1469 FTensor::Tensor2<double, 3, 3> t_d;
1470
1471 double lambda = 1;
1472 if (auto arc_dof = dataAtPts->arcLengthDof.lock()) {
1473 lambda = arc_dof->getFieldData();
1474 }
1475 FTensor1 t_force(dAta.data.data.value3, dAta.data.data.value4,
1476 dAta.data.data.value5);
1477 t_force(i) *= dAta.data.data.value1;
1478 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1479
1480 double nrm2 = std::sqrt(t_normal(i) * t_normal(i));
1481 double a = -0.5 * t_w * nrm2 * lambda;
1482
1483 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1484
1485 int bbc = 0;
1486 for (; bbc != nb_base_fun_col; ++bbc) {
1487
1488 FTensor::Tensor0<double *> t_row_base(&row_data.getN()(gg, 0));
1489
1490 int bbr = 0;
1491 for (; bbr != nb_base_fun_row; ++bbr) {
1492
1493 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
1494
1495 // TODO: handle hoGeometry
1496
1497 t_assemble(i, j) -= a * t_row_base * t_inv_H(l, j) *
1498 t_col_diff_base(m) * t_inv_H(m, i) * t_h(k, l) *
1499 t_force(k);
1500
1501 ++t_row_base;
1502 }
1503 ++t_col_diff_base;
1504 }
1505 ++t_w;
1506 ++t_h;
1507 ++t_inv_H;
1508 ++t_normal;
1509 }
1510
1511 MoFEMFunctionReturn(0);
1512}
1513
1514NeumannForcesSurface::OpNeumannFlux::OpNeumannFlux(
1515 const std::string field_name, Vec _F, bCPressure &data,
1516 boost::ptr_vector<MethodForForceScaling> &methods_op, bool ho_geometry)
1517 : FaceElementForcesAndSourcesCore::UserDataOperator(
1518 field_name, UserDataOperator::OPROW),
1519 F(_F), dAta(data), methodsOp(methods_op), hoGeometry(ho_geometry) {}
1520
1521MoFEMErrorCode
1522NeumannForcesSurface::OpNeumannFlux::doWork(int side, EntityType type,
1523 EntitiesFieldData::EntData &data) {
1524 MoFEMFunctionBegin;
1525
1526 if (data.getIndices().size() == 0)
1527 MoFEMFunctionReturnHot(0);
1528 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1529 dAta.tRis.end()) {
1530 MoFEMFunctionReturnHot(0);
1531 }
1532
1533 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
1534 int nb_row_dofs = data.getIndices().size() / rank;
1535
1536 Nf.resize(data.getIndices().size(), false);
1537 Nf.clear();
1538
1539 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
1540
1541 for (unsigned int gg = 0; gg != data.getN().size1(); ++gg) {
1542
1543 double val = getGaussPts()(2, gg);
1544 double flux;
1545 if (hoGeometry || fe_type == MBQUAD) {
1546 double area = cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
1547 if (fe_type == MBTRI)
1548 area /= 2;
1549 flux = dAta.data.data.value1 * area;
1550 } else {
1551 flux = dAta.data.data.value1 * getArea();
1552 }
1553
1554 cblas_daxpy(nb_row_dofs, val * flux, &data.getN()(gg, 0), 1,
1555 &*Nf.data().begin(), 1);
1556 }
1557
1558 CHKERR MethodForForceScaling::applyScale(getFEMethod(), methodsOp, Nf);
1559
1560 auto get_f = [&]() {
1561 if (F == PETSC_NULL)
1562 return getKSPf();
1563 return F;
1564 };
1565
1566 CHKERR VecSetValues(get_f(), data, &*Nf.data().begin(), ADD_VALUES);
1567
1568 MoFEMFunctionReturn(0);
1569}
1570
1571MoFEMErrorCode NeumannForcesSurface::addForce(const std::string field_name,
1572 Vec F, int ms_id,
1573 bool ho_geometry,
1574 bool block_set) {
1575 const CubitMeshSets *cubit_meshset_ptr;
1576 MeshsetsManager *mmanager_ptr;
1577 MoFEMFunctionBegin;
1578 CHKERR mField.getInterface(mmanager_ptr);
1579 if (block_set) {
1580 // Add data from block set.
1581 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, BLOCKSET,
1582 &cubit_meshset_ptr);
1583 std::vector<double> mydata;
1584 CHKERR cubit_meshset_ptr->getAttributes(mydata);
1585 VectorDouble force(mydata.size());
1586 for (unsigned int ii = 0; ii != mydata.size(); ++ii) {
1587 force[ii] = mydata[ii];
1588 }
1589 // Read forces from BLOCKSET Force (if exists)
1590 if (force.empty()) {
1591 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Force not given");
1592 }
1593 // Assign values from BLOCKSET FORCE to RHS vector. Info about native Cubit
1594 // BC data structure can be found in BCData.hpp
1595 const string name = "Force";
1596 strncpy(mapForce[ms_id].data.data.name, name.c_str(),
1597 name.size() > 5 ? 5 : name.size());
1598 double magnitude = std::sqrt(force[0] * force[0] + force[1] * force[1] +
1599 force[2] * force[2]);
1600 mapForce[ms_id].data.data.value1 = -magnitude; //< Force magnitude
1601 mapForce[ms_id].data.data.value2 = 0;
1602 mapForce[ms_id].data.data.value3 =
1603 force[0] / magnitude; //< X-component of force vector
1604 mapForce[ms_id].data.data.value4 =
1605 force[1] / magnitude; //< Y-component of force vector
1606 mapForce[ms_id].data.data.value5 =
1607 force[2] / magnitude; //< Z-component of force vector
1608 mapForce[ms_id].data.data.value6 = 0;
1609 mapForce[ms_id].data.data.value7 = 0;
1610 mapForce[ms_id].data.data.value8 = 0;
1611 mapForce[ms_id].data.data.zero[0] = 0;
1612 mapForce[ms_id].data.data.zero[1] = 0;
1613 mapForce[ms_id].data.data.zero[2] = 0;
1614 mapForce[ms_id].data.data.zero2 = 0;
1615
1616 CHKERR mField.get_moab().get_entities_by_dimension(
1617 cubit_meshset_ptr->meshset, 2, mapForce[ms_id].tRis, true);
1618 fe.getOpPtrVector().push_back(new OpNeumannForce(
1619 field_name, F, mapForce[ms_id], methodsOp, ho_geometry));
1620
1621 // SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,"Not implemented");
1622 } else {
1623 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, NODESET, &cubit_meshset_ptr);
1624 CHKERR cubit_meshset_ptr->getBcDataStructure(mapForce[ms_id].data);
1625 CHKERR mField.get_moab().get_entities_by_dimension(
1626 cubit_meshset_ptr->meshset, 2, mapForce[ms_id].tRis, true);
1627 fe.getOpPtrVector().push_back(new OpNeumannForce(
1628 field_name, F, mapForce[ms_id], methodsOp, ho_geometry));
1629 }
1630 MoFEMFunctionReturn(0);
1631}
1632
1633MoFEMErrorCode NeumannForcesSurface::addForceAle(
1634 const std::string x_field, const std::string X_field,
1635 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1636 std::string side_fe_name, Vec F, Mat aij, int ms_id, bool ho_geometry,
1637 bool block_set, bool ignore_material_force) {
1638 const CubitMeshSets *cubit_meshset_ptr;
1639 MeshsetsManager *mmanager_ptr;
1640 MoFEMFunctionBegin;
1641 CHKERR mField.getInterface(mmanager_ptr);
1642 if (block_set) {
1643 // Add data from block set.
1644 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, BLOCKSET,
1645 &cubit_meshset_ptr);
1646 std::vector<double> mydata;
1647 CHKERR cubit_meshset_ptr->getAttributes(mydata);
1648 VectorDouble force(mydata.size());
1649 for (unsigned int ii = 0; ii != mydata.size(); ++ii) {
1650 force[ii] = mydata[ii];
1651 }
1652 // Read forces from BLOCKSET Force (if exists)
1653 if (force.empty()) {
1654 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Force not given");
1655 }
1656 // Assign values from BLOCKSET FORCE to RHS vector. Info about native Cubit
1657 // BC data structure can be found in BCData.hpp
1658 const string name = "Force";
1659 strncpy(mapForce[ms_id].data.data.name, name.c_str(),
1660 name.size() > 5 ? 5 : name.size());
1661 double magnitude = std::sqrt(force[0] * force[0] + force[1] * force[1] +
1662 force[2] * force[2]);
1663 mapForce[ms_id].data.data.value1 = -magnitude; //< Force magnitude
1664 mapForce[ms_id].data.data.value2 = 0;
1665 mapForce[ms_id].data.data.value3 =
1666 force[0] / magnitude; //< X-component of force vector
1667 mapForce[ms_id].data.data.value4 =
1668 force[1] / magnitude; //< Y-component of force vector
1669 mapForce[ms_id].data.data.value5 =
1670 force[2] / magnitude; //< Z-component of force vector
1671 mapForce[ms_id].data.data.value6 = 0;
1672 mapForce[ms_id].data.data.value7 = 0;
1673 mapForce[ms_id].data.data.value8 = 0;
1674 mapForce[ms_id].data.data.zero[0] = 0;
1675 mapForce[ms_id].data.data.zero[1] = 0;
1676 mapForce[ms_id].data.data.zero[2] = 0;
1677 mapForce[ms_id].data.data.zero2 = 0;
1678
1679 CHKERR mField.get_moab().get_entities_by_dimension(
1680 cubit_meshset_ptr->meshset, 2, mapForce[ms_id].tRis, true);
1681 // SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,"Not implemented");
1682 } else {
1683 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, NODESET, &cubit_meshset_ptr);
1684 CHKERR cubit_meshset_ptr->getBcDataStructure(mapForce[ms_id].data);
1685 CHKERR mField.get_moab().get_entities_by_dimension(
1686 cubit_meshset_ptr->meshset, 2, mapForce[ms_id].tRis, true);
1687 }
1688
1689 /* LEFT-HAND SIDE (SPATIAL) */
1690 feLhs.getOpPtrVector().push_back(new OpGetTangent(X_field, data_at_pts));
1691
1692 feLhs.getOpPtrVector().push_back(new OpNeumannSurfaceForceLhs_dx_dX(
1693 x_field, X_field, data_at_pts, aij, mapForce[ms_id], ho_geometry));
1694
1695 /* RIGHT-HAND SIDE (MATERIAL) */
1696 if (!ignore_material_force) {
1697 // Side volume element computes the deformation gradient F=hH^-1
1698 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideRhs =
1699 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1700
1701 feMatSideRhs->getOpPtrVector().push_back(
1702 new OpCalculateVectorFieldGradient<3, 3>(X_field,
1703 data_at_pts->getHMatPtr()));
1704 feMatSideRhs->getOpPtrVector().push_back(
1705 new OpCalculateVectorFieldGradient<3, 3>(
1706 x_field, data_at_pts->getSmallhMatPtr()));
1707
1708 feMatSideRhs->getOpPtrVector().push_back(
1709 new OpCalculateDeformation(X_field, data_at_pts, ho_geometry));
1710
1711 feMatRhs.getOpPtrVector().push_back(new OpNeumannSurfaceForceMaterialRhs_dX(
1712 X_field, data_at_pts, feMatSideRhs, side_fe_name, F, mapForce[ms_id],
1713 ho_geometry));
1714
1715 /* LEFT-HAND SIDE (MATERIAL) */
1716
1717 // Side volume element computes linearisation with spatial coordinates
1718 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dx =
1719 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1720 // Side volume element computes linearisation with material coordinates
1721 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dX =
1722 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1723
1724 feMatSideLhs_dx->getOpPtrVector().push_back(
1725 new OpCalculateVectorFieldGradient<3, 3>(X_field,
1726 data_at_pts->getHMatPtr()));
1727 feMatSideLhs_dx->getOpPtrVector().push_back(
1728 new OpCalculateVectorFieldGradient<3, 3>(
1729 x_field, data_at_pts->getSmallhMatPtr()));
1730
1731 feMatSideLhs_dx->getOpPtrVector().push_back(
1732 new OpCalculateDeformation(X_field, data_at_pts, ho_geometry));
1733 feMatSideLhs_dx->getOpPtrVector().push_back(
1734 new OpNeumannSurfaceForceMaterialVolOnSideLhs_dX_dx(
1735 X_field, x_field, data_at_pts, aij, mapForce[ms_id], ho_geometry));
1736
1737 feMatLhs.getOpPtrVector().push_back(
1738 new OpNeumannSurfaceForceMaterialLhs_dX_dx(
1739 X_field, x_field, data_at_pts, feMatSideLhs_dx, side_fe_name, aij,
1740 mapForce[ms_id], ho_geometry));
1741
1742 feMatSideLhs_dX->getOpPtrVector().push_back(
1743 new OpCalculateVectorFieldGradient<3, 3>(X_field,
1744 data_at_pts->getHMatPtr()));
1745 feMatSideLhs_dX->getOpPtrVector().push_back(
1746 new OpCalculateVectorFieldGradient<3, 3>(
1747 x_field, data_at_pts->getSmallhMatPtr()));
1748 feMatSideLhs_dX->getOpPtrVector().push_back(
1749 new OpCalculateDeformation(X_field, data_at_pts, ho_geometry));
1750 feMatSideLhs_dX->getOpPtrVector().push_back(
1751 new OpNeumannSurfaceForceMaterialVolOnSideLhs_dX_dX(
1752 X_field, X_field, data_at_pts, aij, mapForce[ms_id], ho_geometry));
1753
1754 feMatLhs.getOpPtrVector().push_back(new OpGetTangent(X_field, data_at_pts));
1755 feMatLhs.getOpPtrVector().push_back(
1756 new OpNeumannSurfaceForceMaterialLhs_dX_dX(
1757 X_field, X_field, data_at_pts, feMatSideLhs_dX, side_fe_name, aij,
1758 mapForce[ms_id], ho_geometry));
1759 }
1760
1761 MoFEMFunctionReturn(0);
1762}
1763
1764MoFEMErrorCode NeumannForcesSurface::addPressure(const std::string field_name,
1765 Vec F, int ms_id,
1766 bool ho_geometry,
1767 bool block_set) {
1768
1769 const CubitMeshSets *cubit_meshset_ptr;
1770 MeshsetsManager *mmanager_ptr;
1771 MoFEMFunctionBegin;
1772 CHKERR mField.getInterface(mmanager_ptr);
1773 if (block_set) {
1774 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, BLOCKSET,
1775 &cubit_meshset_ptr);
1776 std::vector<double> mydata;
1777 CHKERR cubit_meshset_ptr->getAttributes(mydata);
1778 VectorDouble pressure(mydata.size());
1779 for (unsigned int ii = 0; ii != mydata.size(); ++ii) {
1780 pressure[ii] = mydata[ii];
1781 }
1782 if (pressure.empty()) {
1783 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Pressure not given");
1784 }
1785 const string name = "Pressure";
1786 strncpy(mapPressure[ms_id].data.data.name, name.c_str(),
1787 name.size() > 8 ? 8 : name.size());
1788 mapPressure[ms_id].data.data.flag1 = 0;
1789 mapPressure[ms_id].data.data.flag2 = 1;
1790 mapPressure[ms_id].data.data.value1 = pressure[0];
1791 mapPressure[ms_id].data.data.zero = 0;
1792 CHKERR mField.get_moab().get_entities_by_dimension(
1793 cubit_meshset_ptr->meshset, 2, mapPressure[ms_id].tRis, true);
1794 fe.getOpPtrVector().push_back(new OpNeumannPressure(
1795 field_name, F, mapPressure[ms_id], methodsOp, ho_geometry));
1796 } else {
1797 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, SIDESET, &cubit_meshset_ptr);
1798 CHKERR cubit_meshset_ptr->getBcDataStructure(mapPressure[ms_id].data);
1799 CHKERR mField.get_moab().get_entities_by_dimension(
1800 cubit_meshset_ptr->meshset, 2, mapPressure[ms_id].tRis, true);
1801 fe.getOpPtrVector().push_back(new OpNeumannPressure(
1802 field_name, F, mapPressure[ms_id], methodsOp, ho_geometry));
1803 }
1804 MoFEMFunctionReturn(0);
1805}
1806
1807MoFEMErrorCode NeumannForcesSurface::addPressureAle(
1808 const std::string x_field, const std::string X_field,
1809 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1810 std::string side_fe_name, Vec F, Mat aij, int ms_id, bool ho_geometry,
1811 bool block_set) {
1812
1813 const CubitMeshSets *cubit_meshset_ptr;
1814 MeshsetsManager *mmanager_ptr;
1815 MoFEMFunctionBegin;
1816 CHKERR mField.getInterface(mmanager_ptr);
1817 if (block_set) {
1818 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, BLOCKSET,
1819 &cubit_meshset_ptr);
1820 std::vector<double> mydata;
1821 CHKERR cubit_meshset_ptr->getAttributes(mydata);
1822 VectorDouble pressure(mydata.size());
1823 for (unsigned int ii = 0; ii != mydata.size(); ++ii) {
1824 pressure[ii] = mydata[ii];
1825 }
1826 if (pressure.empty()) {
1827 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Pressure not given");
1828 }
1829 const string name = "Pressure";
1830 strncpy(mapPressure[ms_id].data.data.name, name.c_str(),
1831 name.size() > 8 ? 8 : name.size());
1832 mapPressure[ms_id].data.data.flag1 = 0;
1833 mapPressure[ms_id].data.data.flag2 = 1;
1834 mapPressure[ms_id].data.data.value1 = pressure[0];
1835 mapPressure[ms_id].data.data.zero = 0;
1836 CHKERR mField.get_moab().get_entities_by_dimension(
1837 cubit_meshset_ptr->meshset, 2, mapPressure[ms_id].tRis, true);
1838
1839 } else {
1840 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, SIDESET, &cubit_meshset_ptr);
1841 CHKERR cubit_meshset_ptr->getBcDataStructure(mapPressure[ms_id].data);
1842 CHKERR mField.get_moab().get_entities_by_dimension(
1843 cubit_meshset_ptr->meshset, 2, mapPressure[ms_id].tRis, true);
1844 }
1845
1846 /* LEFT-HAND SIDE (SPATIAL) */
1847
1848 feLhs.getOpPtrVector().push_back(new OpGetTangent(X_field, data_at_pts));
1849
1850 feLhs.getOpPtrVector().push_back(new OpNeumannPressureLhs_dx_dX(
1851 x_field, X_field, data_at_pts, aij, mapPressure[ms_id], ho_geometry));
1852
1853 /* RIGHT-HAND SIDE (MATERIAL) */
1854
1855 // Side volume element computes the deformation gradient F=hH^-1
1856 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideRhs =
1857 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1858
1859 feMatSideRhs->getOpPtrVector().push_back(
1860 new OpCalculateVectorFieldGradient<3, 3>(X_field,
1861 data_at_pts->getHMatPtr()));
1862 feMatSideRhs->getOpPtrVector().push_back(
1863 new OpCalculateVectorFieldGradient<3, 3>(x_field,
1864 data_at_pts->getSmallhMatPtr()));
1865
1866 feMatSideRhs->getOpPtrVector().push_back(
1867 new OpCalculateDeformation(X_field, data_at_pts, ho_geometry));
1868
1869 feMatRhs.getOpPtrVector().push_back(new OpNeumannPressureMaterialRhs_dX(
1870 X_field, data_at_pts, feMatSideRhs, side_fe_name, F, mapPressure[ms_id],
1871 ho_geometry));
1872
1873 /* LEFT-HAND SIDE (MATERIAL) */
1874
1875 // Side volume element computes linearisation with spatial coordinates
1876 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dx =
1877 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1878 // Side volume element computes linearisation with material coordinates
1879 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dX =
1880 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1881
1882 feMatSideLhs_dx->getOpPtrVector().push_back(
1883 new OpCalculateVectorFieldGradient<3, 3>(X_field,
1884 data_at_pts->getHMatPtr()));
1885 feMatSideLhs_dx->getOpPtrVector().push_back(
1886 new OpCalculateVectorFieldGradient<3, 3>(x_field,
1887 data_at_pts->getSmallhMatPtr()));
1888
1889 feMatSideLhs_dx->getOpPtrVector().push_back(
1890 new OpCalculateDeformation(X_field, data_at_pts, ho_geometry));
1891 feMatSideLhs_dx->getOpPtrVector().push_back(
1892 new OpNeumannPressureMaterialVolOnSideLhs_dX_dx(
1893 X_field, x_field, data_at_pts, aij, mapPressure[ms_id], ho_geometry));
1894
1895 feMatLhs.getOpPtrVector().push_back(new OpNeumannPressureMaterialLhs_dX_dx(
1896 X_field, x_field, data_at_pts, feMatSideLhs_dx, side_fe_name, aij,
1897 mapPressure[ms_id], ho_geometry));
1898
1899 feMatSideLhs_dX->getOpPtrVector().push_back(
1900 new OpCalculateVectorFieldGradient<3, 3>(X_field,
1901 data_at_pts->getHMatPtr()));
1902 feMatSideLhs_dX->getOpPtrVector().push_back(
1903 new OpCalculateVectorFieldGradient<3, 3>(x_field,
1904 data_at_pts->getSmallhMatPtr()));
1905 feMatSideLhs_dX->getOpPtrVector().push_back(
1906 new OpCalculateDeformation(X_field, data_at_pts, ho_geometry));
1907 feMatSideLhs_dX->getOpPtrVector().push_back(
1908 new OpNeumannPressureMaterialVolOnSideLhs_dX_dX(
1909 X_field, X_field, data_at_pts, aij, mapPressure[ms_id], ho_geometry));
1910
1911 feMatLhs.getOpPtrVector().push_back(new OpGetTangent(X_field, data_at_pts));
1912 feMatLhs.getOpPtrVector().push_back(new OpNeumannPressureMaterialLhs_dX_dX(
1913 X_field, X_field, data_at_pts, feMatSideLhs_dX, side_fe_name, aij,
1914 mapPressure[ms_id], ho_geometry));
1915
1916 MoFEMFunctionReturn(0);
1917}
1918
1919MoFEMErrorCode
1920NeumannForcesSurface::addLinearPressure(const std::string field_name, Vec F,
1921 int ms_id, bool ho_geometry) {
1922
1923 const CubitMeshSets *cubit_meshset_ptr;
1924 MeshsetsManager *mmanager_ptr;
1925 MoFEMFunctionBegin;
1926 CHKERR mField.getInterface(mmanager_ptr);
1927 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, BLOCKSET, &cubit_meshset_ptr);
1928 std::vector<double> mydata;
1929 CHKERR cubit_meshset_ptr->getAttributes(mydata);
1930 if (mydata.size() < 4)
1931 SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1932 "Should be four block attributes but is %d", mydata.size());
1933 VectorDouble3 pressure_coeffs(3);
1934 for (unsigned int ii = 0; ii != 3; ++ii) {
1935 pressure_coeffs[ii] = mydata[ii];
1936 }
1937 const double pressure_shift = mydata[3];
1938
1939 Range tris;
1940 CHKERR mField.get_moab().get_entities_by_dimension(cubit_meshset_ptr->meshset,
1941 2, tris, true);
1942 boost::shared_ptr<MethodForAnalyticalForce> analytical_force_op(
1943 new LinearVaringPresssure(pressure_coeffs, pressure_shift));
1944 fe.getOpPtrVector().push_back(new OpNeumannForceAnalytical(
1945 field_name, F, tris, methodsOp, analytical_force_op, ho_geometry));
1946
1947 MoFEMFunctionReturn(0);
1948}
1949
1950MoFEMErrorCode NeumannForcesSurface::addPreassure(const std::string field_name,
1951 Vec F, int ms_id,
1952 bool ho_geometry,
1953 bool block_set) {
1954 return NeumannForcesSurface::addPressure(field_name, F, ms_id, ho_geometry,
1955 block_set);
1956}
1957
1958MoFEMErrorCode NeumannForcesSurface::addFlux(const std::string field_name,
1959 Vec F, int ms_id,
1960 bool ho_geometry) {
1961 const CubitMeshSets *cubit_meshset_ptr;
1962 MeshsetsManager *mmanager_ptr;
1963 MoFEMFunctionBegin;
1964 CHKERR mField.getInterface(mmanager_ptr);
1965 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, SIDESET, &cubit_meshset_ptr);
1966 CHKERR cubit_meshset_ptr->getBcDataStructure(mapPressure[ms_id].data);
1967 CHKERR mField.get_moab().get_entities_by_dimension(
1968 cubit_meshset_ptr->meshset, 2, mapPressure[ms_id].tRis, true);
1969 fe.getOpPtrVector().push_back(new OpNeumannFlux(
1970 field_name, F, mapPressure[ms_id], methodsOp, ho_geometry));
1971 MoFEMFunctionReturn(0);
1972}
1973
1974MoFEMErrorCode MetaNeumannForces::addNeumannBCElements(
1975 MoFEM::Interface &m_field, const std::string field_name,
1976 const std::string mesh_nodals_positions, Range *intersect_ptr) {
1977 MoFEMFunctionBegin;
1978
1979 // Define boundary element that operates on rows, columns and data of a
1980 // given field
1981 CHKERR m_field.add_finite_element("FORCE_FE", MF_ZERO);
1982 CHKERR m_field.modify_finite_element_add_field_row("FORCE_FE", field_name);
1983 CHKERR m_field.modify_finite_element_add_field_col("FORCE_FE", field_name);
1984 CHKERR m_field.modify_finite_element_add_field_data("FORCE_FE", field_name);
1985 if (m_field.check_field(mesh_nodals_positions)) {
1986 CHKERR m_field.modify_finite_element_add_field_data("FORCE_FE",
1987 mesh_nodals_positions);
1988 }
1989 // Add entities to that element, here we add all triangles with FORCESET
1990 // from cubit
1991 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
1992 it)) {
1993 Range tris;
1994 CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 2, tris,
1995 true);
1996 if (intersect_ptr)
1997 tris = intersect(tris, *intersect_ptr);
1998 CHKERR m_field.add_ents_to_finite_element_by_dim(tris, 2, "FORCE_FE");
1999 }
2000
2001 CHKERR m_field.add_finite_element("PRESSURE_FE", MF_ZERO);
2002 CHKERR m_field.modify_finite_element_add_field_row("PRESSURE_FE", field_name);
2003 CHKERR m_field.modify_finite_element_add_field_col("PRESSURE_FE", field_name);
2004 CHKERR m_field.modify_finite_element_add_field_data("PRESSURE_FE",
2005 field_name);
2006 if (m_field.check_field(mesh_nodals_positions)) {
2007 CHKERR m_field.modify_finite_element_add_field_data("PRESSURE_FE",
2008 mesh_nodals_positions);
2009 }
2010
2011 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,
2012 SIDESET | PRESSURESET, it)) {
2013 Range tris;
2014 CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 2, tris,
2015 true);
2016 if (intersect_ptr)
2017 tris = intersect(tris, *intersect_ptr);
2018 CHKERR m_field.add_ents_to_finite_element_by_dim(tris, 2, "PRESSURE_FE");
2019 }
2020
2021 // Reading forces from BLOCKSET
2022
2023 const string block_set_force_name("FORCE");
2024 // search for block named FORCE and add its attributes to FORCE_FE element
2025 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
2026 if (it->getName().compare(0, block_set_force_name.length(),
2027 block_set_force_name) == 0) {
2028 Range tris;
2029 CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 2, tris,
2030 true);
2031 if (intersect_ptr)
2032 tris = intersect(tris, *intersect_ptr);
2033 CHKERR m_field.add_ents_to_finite_element_by_dim(tris, 2, "FORCE_FE");
2034 }
2035 }
2036 // search for block named PRESSURE and add its attributes to PRESSURE_FE
2037 // element
2038 const string block_set_pressure_name("PRESSURE");
2039 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
2040 if (it->getName().compare(0, block_set_pressure_name.length(),
2041 block_set_pressure_name) == 0) {
2042 Range tris;
2043 CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 2, tris,
2044 true);
2045 if (intersect_ptr)
2046 tris = intersect(tris, *intersect_ptr);
2047 CHKERR m_field.add_ents_to_finite_element_by_dim(tris, 2, "PRESSURE_FE");
2048 }
2049 }
2050
2051 // search for block named LINEAR_PRESSURE and add its attributes to
2052 // PRESSURE_FE element
2053 const string block_set_linear_pressure_name("LINEAR_PRESSURE");
2054 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
2055 if (it->getName().compare(0, block_set_linear_pressure_name.length(),
2056 block_set_linear_pressure_name) == 0) {
2057 Range tris;
2058 CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 2, tris,
2059 true);
2060 if (intersect_ptr)
2061 tris = intersect(tris, *intersect_ptr);
2062 CHKERR m_field.add_ents_to_finite_element_by_dim(tris, 2, "PRESSURE_FE");
2063 }
2064 }
2065
2066 MoFEMFunctionReturn(0);
2067}
2068
2069MoFEMErrorCode MetaNeumannForces::setMomentumFluxOperators(
2070 MoFEM::Interface &m_field,
2071 boost::ptr_map<std::string, NeumannForcesSurface> &neumann_forces, Vec F,
2072 const std::string field_name, const std::string mesh_nodals_positions) {
2073 MoFEMFunctionBegin;
2074 bool ho_geometry = m_field.check_field(mesh_nodals_positions);
2075 // Add forces
2076 {
2077 std::string fe_name = "FORCE_FE";
2078 neumann_forces.insert(fe_name, new NeumannForcesSurface(m_field));
2079 auto &nf = neumann_forces.at(fe_name);
2080 auto &fe = nf.getLoopFe();
2081
2082 if (ho_geometry)
2083 CHKERR addHOOpsFace3D(mesh_nodals_positions, fe, false, false);
2084
2085 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
2086 it)) {
2087 CHKERR nf.addForce(field_name, F, it->getMeshsetId(), ho_geometry, false);
2088 }
2089 // Reading forces from BLOCKSET
2090 const string block_set_force_name("FORCE");
2091 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
2092 if (it->getName().compare(0, block_set_force_name.length(),
2093 block_set_force_name) == 0) {
2094 CHKERR nf.addForce(field_name, F, it->getMeshsetId(), ho_geometry,
2095 true);
2096 }
2097 }
2098 }
2099
2100 // Add pressures
2101 {
2102 std::string fe_name = "PRESSURE_FE";
2103 neumann_forces.insert(fe_name, new NeumannForcesSurface(m_field));
2104
2105 auto &nf = neumann_forces.at(fe_name);
2106 auto &fe = nf.getLoopFe();
2107
2108 if (ho_geometry)
2109 CHKERR addHOOpsFace3D(mesh_nodals_positions, fe, false, false);
2110
2111 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
2112 m_field, SIDESET | PRESSURESET, it)) {
2113 CHKERR nf.addPressure(field_name, F, it->getMeshsetId(), ho_geometry,
2114 false);
2115 }
2116
2117 // Reading pressures from BLOCKSET
2118 const string block_set_pressure_name("PRESSURE");
2119 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
2120 if (it->getName().compare(0, block_set_pressure_name.length(),
2121 block_set_pressure_name) == 0) {
2122 CHKERR nf.addPressure(field_name, F, it->getMeshsetId(), ho_geometry,
2123 true);
2124 }
2125 }
2126
2127 // Reading pressures from BLOCKSET
2128 const string block_set_linear_pressure_name("LINEAR_PRESSURE");
2129 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
2130 if (it->getName().compare(0, block_set_linear_pressure_name.length(),
2131 block_set_linear_pressure_name) == 0) {
2132 CHKERR nf.addLinearPressure(field_name, F, it->getMeshsetId(),
2133 ho_geometry);
2134 }
2135 }
2136 }
2137
2138 MoFEMFunctionReturn(0);
2139}
2140
2141MoFEMErrorCode MetaNeumannForces::addNeumannFluxBCElements(
2142 MoFEM::Interface &m_field, const std::string field_name,
2143 const std::string mesh_nodals_positions) {
2144 MoFEMFunctionBegin;
2145
2146 CHKERR m_field.add_finite_element("FLUX_FE", MF_ZERO);
2147 CHKERR m_field.modify_finite_element_add_field_row("FLUX_FE", field_name);
2148 CHKERR m_field.modify_finite_element_add_field_col("FLUX_FE", field_name);
2149 CHKERR m_field.modify_finite_element_add_field_data("FLUX_FE", field_name);
2150 if (m_field.check_field(mesh_nodals_positions)) {
2151 CHKERR m_field.modify_finite_element_add_field_data("FLUX_FE",
2152 mesh_nodals_positions);
2153 }
2154
2155 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,
2156 SIDESET | PRESSURESET, it)) {
2157 Range tris;
2158 CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 2, tris,
2159 true);
2160 CHKERR m_field.add_ents_to_finite_element_by_dim(tris, 2, "FLUX_FE");
2161 }
2162
2163 MoFEMFunctionReturn(0);
2164}
2165
2166MoFEMErrorCode MetaNeumannForces::setMassFluxOperators(
2167 MoFEM::Interface &m_field,
2168 boost::ptr_map<std::string, NeumannForcesSurface> &neumann_forces, Vec F,
2169 const std::string field_name, const std::string mesh_nodals_positions) {
2170 MoFEMFunctionBegin;
2171 bool ho_geometry = m_field.check_field(mesh_nodals_positions);
2172
2173 string fe_name;
2174 fe_name = "FLUX_FE";
2175 neumann_forces.insert(fe_name, new NeumannForcesSurface(m_field));
2176
2177 auto &nf = neumann_forces.at(fe_name);
2178 auto &fe = nf.getLoopFe();
2179
2180 if (ho_geometry)
2181 CHKERR addHOOpsFace3D(mesh_nodals_positions, fe, false, false);
2182
2183 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,
2184 SIDESET | PRESSURESET, it)) {
2185 CHKERR nf.addFlux(field_name, F, it->getMeshsetId(), ho_geometry);
2186 }
2187 MoFEMFunctionReturn(0);
2188}
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
a
constexpr double a
Definition: approx_sphere.cpp:30
FTensor::Tensor1::data
T data[Tensor_Dim]
Definition: Tensor1_value.hpp:10
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
PRESSURESET
@ PRESSURESET
Definition: definitions.h:152
FORCESET
@ FORCESET
Definition: definitions.h:151
NODESET
@ NODESET
Definition: definitions.h:146
SIDESET
@ SIDESET
Definition: definitions.h:147
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
m
FTensor::Index< 'm', SPACE_DIM > m
Definition: fluid_structure_eigenproblem.cpp:65
F
@ F
Definition: free_surface.cpp:394
get_f
auto get_f
Definition: free_surface.cpp:221
MoFEM::CoreInterface::add_ents_to_finite_element_by_dim
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
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:563
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::VectorDouble3
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::addHOOpsFace3D
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Definition: HODataOperators.hpp:779
MoFEM::invertTensor3by3
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1608
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:1631
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1557
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
MoFEM::VectorDofs
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
Definition: EntitiesFieldData.hpp:23
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
_F
#define _F(n)
Definition: quad.c:25
MetaNeumannForces::setMassFluxOperators
static MoFEMErrorCode setMassFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Definition: SurfacePressure.cpp:2166
MetaNeumannForces::addNeumannFluxBCElements
static MoFEMErrorCode addNeumannFluxBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Definition: SurfacePressure.cpp:2141
MetaNeumannForces::addNeumannBCElements
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
Definition: SurfacePressure.cpp:1974
MetaNeumannForces::setMomentumFluxOperators
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
Definition: SurfacePressure.cpp:2069
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::CubitMeshSets
this struct keeps basic methods for moab meshset about material and boundary conditions
Definition: BCMultiIndices.hpp:19
MoFEM::CubitMeshSets::getAttributes
MoFEMErrorCode getAttributes(std::vector< double > &attributes) const
get Cubit block attributes
Definition: BCMultiIndices.cpp:274
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:127
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:387
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1489
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:1305
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1241
MoFEM::EntitiesFieldData::EntData::getFTensor1FieldData
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
Definition: EntitiesFieldData.hpp:1275
MoFEM::EntitiesFieldData::EntData::getFieldDofs
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
Definition: EntitiesFieldData.hpp:1256
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1201
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::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1558
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
NeumannForcesSurface::LinearVaringPresssure::getForce
MoFEMErrorCode getForce(const EntityHandle ent, const VectorDouble3 &coords, const VectorDouble3 &normal, VectorDouble3 &force)
Definition: SurfacePressure.cpp:14
NeumannForcesSurface::MyTriangleFE::MyTriangleFE
MyTriangleFE(MoFEM::Interface &m_field)
Definition: SurfacePressure.cpp:23
NeumannForcesSurface::OpCalculateDeformation
Operator for computing deformation gradients in side volumes.
Definition: SurfacePressure.hpp:401
NeumannForcesSurface::OpCalculateDeformation::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfacePressure.cpp:470
NeumannForcesSurface::OpGetTangent
Operator for computing tangent vectors.
Definition: SurfacePressure.hpp:239
NeumannForcesSurface::OpGetTangent::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfacePressure.cpp:226
NeumannForcesSurface::OpNeumannFlux::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfacePressure.cpp:1522
NeumannForcesSurface::OpNeumannFlux::OpNeumannFlux
OpNeumannFlux(const std::string field_name, Vec _F, bCPressure &data, boost::ptr_vector< MethodForForceScaling > &methods_op, bool ho_geometry)
Definition: SurfacePressure.cpp:1514
NeumannForcesSurface::OpNeumannForceAnalytical::OpNeumannForceAnalytical
OpNeumannForceAnalytical(const std::string field_name, Vec f, const Range tris, boost::ptr_vector< MethodForForceScaling > &methods_op, boost::shared_ptr< MethodForAnalyticalForce > &analytical_force_op, const bool ho_geometry=false)
Definition: SurfacePressure.cpp:97
NeumannForcesSurface::OpNeumannForceAnalytical::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfacePressure.cpp:107
NeumannForcesSurface::OpNeumannForce::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate surface force (traction)
Definition: SurfacePressure.cpp:34
NeumannForcesSurface::OpNeumannForce::OpNeumannForce
OpNeumannForce(const std::string field_name, Vec _F, bCForce &data, boost::ptr_vector< MethodForForceScaling > &methods_op, bool ho_geometry=false)
Definition: SurfacePressure.cpp:26
NeumannForcesSurface::OpNeumannPressure
RHS-operator for pressure element (spatial configuration)
Definition: SurfacePressure.hpp:168
NeumannForcesSurface::OpNeumannPressure::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate pressure.
Definition: SurfacePressure.cpp:174
NeumannForcesSurface::OpNeumannPressure::OpNeumannPressure
OpNeumannPressure(const std::string field_name, Vec _F, bCPressure &data, boost::ptr_vector< MethodForForceScaling > &methods_op, bool ho_geometry=false)
Definition: SurfacePressure.cpp:167
NeumannForcesSurface::OpNeumannPressureLhs_dx_dX
LHS-operator for pressure element (spatial configuration)
Definition: SurfacePressure.hpp:259
NeumannForcesSurface::OpNeumannPressureLhs_dx_dX::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Compute left-hand side.
Definition: SurfacePressure.cpp:265
NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dX
LHS-operator for the pressure element (material configuration)
Definition: SurfacePressure.hpp:717
NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dX::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: SurfacePressure.cpp:750
NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Compute part of the left-hand side.
Definition: SurfacePressure.cpp:793
NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dx
LHS-operator for the pressure element (material configuration)
Definition: SurfacePressure.hpp:831
NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dx::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: SurfacePressure.cpp:1063
NeumannForcesSurface::OpNeumannPressureMaterialLhs::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Definition: SurfacePressure.cpp:1030
NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX
RHS-operator for the pressure element (material configuration)
Definition: SurfacePressure.hpp:425
NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data)
Definition: SurfacePressure.cpp:537
NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data)
Definition: SurfacePressure.cpp:583
NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &row_data)
Integrate pressure in the material configuration.
Definition: SurfacePressure.cpp:506
NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs_dX_dX
LHS-operator (material configuration) on the side volume.
Definition: SurfacePressure.hpp:1068
NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs_dX_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
Definition: SurfacePressure.cpp:1376
NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs_dX_dx
LHS-operator (material configuration) on the side volume.
Definition: SurfacePressure.hpp:989
NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs_dX_dx::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
Definition: SurfacePressure.cpp:1183
NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: SurfacePressure.cpp:1107
NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Definition: SurfacePressure.cpp:1308
NeumannForcesSurface::OpNeumannSurfaceForceLhs_dx_dX
LHS-operator for surface force element (spatial configuration)
Definition: SurfacePressure.hpp:326
NeumannForcesSurface::OpNeumannSurfaceForceLhs_dx_dX::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Compute left-hand side.
Definition: SurfacePressure.cpp:363
NeumannForcesSurface::OpNeumannSurfaceForceMaterialLhs_dX_dX
LHS-operator for the surface force element (material configuration)
Definition: SurfacePressure.hpp:774
NeumannForcesSurface::OpNeumannSurfaceForceMaterialLhs_dX_dX::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: SurfacePressure.cpp:903
NeumannForcesSurface::OpNeumannSurfaceForceMaterialLhs_dX_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Compute part of the left-hand side.
Definition: SurfacePressure.cpp:946
NeumannForcesSurface::OpNeumannSurfaceForceMaterialLhs_dX_dx
LHS-operator for the surface force element (material configuration)
Definition: SurfacePressure.hpp:862
NeumannForcesSurface::OpNeumannSurfaceForceMaterialLhs_dX_dx::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: SurfacePressure.cpp:1085
NeumannForcesSurface::OpNeumannSurfaceForceMaterialLhs::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Definition: SurfacePressure.cpp:869
NeumannForcesSurface::OpNeumannSurfaceForceMaterialRhs_dX
RHS-operator for the surface force element (material configuration)
Definition: SurfacePressure.hpp:496
NeumannForcesSurface::OpNeumannSurfaceForceMaterialRhs_dX::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data)
Definition: SurfacePressure.cpp:709
NeumannForcesSurface::OpNeumannSurfaceForceMaterialRhs_dX::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &row_data)
Integrate surface force in the material configuration.
Definition: SurfacePressure.cpp:625
NeumannForcesSurface::OpNeumannSurfaceForceMaterialRhs_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data)
Definition: SurfacePressure.cpp:657
NeumannForcesSurface::OpNeumannSurfaceForceMaterialVolOnSideLhs_dX_dX
LHS-operator (material configuration) on the side volume.
Definition: SurfacePressure.hpp:1108
NeumannForcesSurface::OpNeumannSurfaceForceMaterialVolOnSideLhs_dX_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
Definition: SurfacePressure.cpp:1445
NeumannForcesSurface::OpNeumannSurfaceForceMaterialVolOnSideLhs_dX_dx
LHS-operator (material configuration) on the side volume.
Definition: SurfacePressure.hpp:1028
NeumannForcesSurface::OpNeumannSurfaceForceMaterialVolOnSideLhs_dX_dx::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
Definition: SurfacePressure.cpp:1244
NeumannForcesSurface::OpNeumannSurfaceForceMaterialVolOnSideLhs::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Definition: SurfacePressure.cpp:1341
NeumannForcesSurface::OpNeumannSurfaceForceMaterialVolOnSideLhs::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: SurfacePressure.cpp:1145
NeumannForcesSurface
Finite element and operators to apply force/pressures applied to surfaces.
Definition: SurfacePressure.hpp:14
NeumannForcesSurface::mapPressure
std::map< int, bCPressure > mapPressure
Definition: SurfacePressure.hpp:95
NeumannForcesSurface::addPressureAle
MoFEMErrorCode addPressureAle(const std::string field_name_1, const std::string field_name_2, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, std::string side_fe_name, Vec F, Mat aij, int ms_id, bool ho_geometry=false, bool block_set=false)
Add operator to calculate pressure on element (in ALE)
Definition: SurfacePressure.cpp:1807
NeumannForcesSurface::addForceAle
MoFEMErrorCode addForceAle(const std::string field_name_1, const std::string field_name_2, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, std::string side_fe_name, Vec F, Mat aij, int ms_id, bool ho_geometry=false, bool block_set=false, bool ignore_material_force=false)
Add operator to calculate forces on element (in ALE)
Definition: SurfacePressure.cpp:1633
NeumannForcesSurface::mapForce
std::map< int, bCForce > mapForce
Definition: SurfacePressure.hpp:89
NeumannForcesSurface::addFlux
MoFEMErrorCode addFlux(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false)
Add flux element operator (integration on face)
Definition: SurfacePressure.cpp:1958
NeumannForcesSurface::methodsOp
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: SurfacePressure.hpp:97
NeumannForcesSurface::addPreassure
DEPRECATED MoFEMErrorCode addPreassure(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false, bool block_set=false)
Definition: SurfacePressure.cpp:1950
NeumannForcesSurface::addForce
MoFEMErrorCode addForce(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false, bool block_set=false)
Add operator to calculate forces on element.
Definition: SurfacePressure.cpp:1571
NeumannForcesSurface::addLinearPressure
MoFEMErrorCode addLinearPressure(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false)
Add operator to calculate pressure on element.
Definition: SurfacePressure.cpp:1920
NeumannForcesSurface::addPressure
MoFEMErrorCode addPressure(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false, bool block_set=false)
Add operator to calculate pressure on element.
Definition: SurfacePressure.cpp:1764
NeumannForcesSurface::analyticalForceOp
boost::ptr_vector< MethodForAnalyticalForce > analyticalForceOp
Definition: SurfacePressure.hpp:98
NeumannForcesSurface::mField
MoFEM::Interface & mField
Definition: SurfacePressure.hpp:16
Generated by Doxygen 1.9.5 and hosted at University of Glasgow