v0.13.1
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/* This file is part of MoFEM.
9 * MoFEM is free software: you can redistribute it and/or modify it under
10 * the terms of the GNU Lesser General Public License as published by the
11 * Free Software Foundation, either version 3 of the License, or (at your
12 * option) any later version.
13 *
14 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
15 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 * License for more details.
18 *
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
21
22using namespace MoFEM;
23
24using namespace boost::numeric;
25
27 const EntityHandle ent, const VectorDouble3 &coords,
28 const VectorDouble3 &normal, VectorDouble3 &force) {
30 const double p = inner_prod(coords, linearConstants) + pressureShift;
31 force = normal * p / norm_2(normal);
33}
34
36 : FaceElementForcesAndSourcesCore(m_field), addToRule(1) {}
37
39 const std::string field_name, Vec _F, bCForce &data,
40 boost::ptr_vector<MethodForForceScaling> &methods_op, bool ho_geometry)
43 F(_F), dAta(data), methodsOp(methods_op), hoGeometry(ho_geometry) {}
44
48
49 if (data.getIndices().size() == 0)
51 EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
52 if (dAta.tRis.find(ent) == dAta.tRis.end())
54
55 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
56 int nb_row_dofs = data.getIndices().size() / rank;
57
58 Nf.resize(data.getIndices().size(), false);
59 Nf.clear();
60
61 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
62
63 for (unsigned int gg = 0; gg != data.getN().size1(); ++gg) {
64
65 // get integration weight and Jacobian of integration point (area of face)
66 double val = getGaussPts()(2, gg);
67 if (hoGeometry || fe_type == MBQUAD) {
68 val *= cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
69 if (fe_type == MBTRI)
70 val /= 2;
71 } else
72 val *= getArea();
73
74 // use data from module
75 for (int rr = 0; rr != rank; ++rr) {
76 double force;
77 if (rr == 0)
78 force = dAta.data.data.value3;
79 else if (rr == 1)
80 force = dAta.data.data.value4;
81 else if (rr == 2)
82 force = dAta.data.data.value5;
83 else
84 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
85 "data inconsistency");
86
87 force *= dAta.data.data.value1;
88 cblas_daxpy(nb_row_dofs, val * force, &data.getN()(gg, 0), 1, &Nf[rr],
89 rank);
90 }
91 }
92
93 // Scale force using user defined scaling operator
95
96 auto get_f = [&]() {
97 if (F == PETSC_NULL)
98 return getKSPf();
99 return F;
100 };
101
102 // Assemble force into vector
103 CHKERR VecSetValues(get_f(), data, &*Nf.data().begin(), ADD_VALUES);
104
106}
107
109 const std::string field_name, Vec f, const Range tris,
110 boost::ptr_vector<MethodForForceScaling> &methods_op,
111 boost::shared_ptr<MethodForAnalyticalForce> &analytical_force_op,
112 const bool ho_geometry)
115 F(f), tRis(tris), methodsOp(methods_op),
116 analyticalForceOp(analytical_force_op), hoGeometry(ho_geometry) {}
117
119 int side, EntityType type, EntitiesFieldData::EntData &data) {
121
122 if (data.getIndices().size() == 0)
124 const EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
125 if (tRis.find(ent) == tRis.end())
127
128 const int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
129 const int nb_row_dofs = data.getIndices().size() / rank;
130
131 nF.resize(data.getIndices().size(), false);
132 nF.clear();
133
134 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
135
136 VectorDouble3 coords(3);
138 VectorDouble3 force(3);
139
140 for (unsigned int gg = 0; gg != data.getN().size1(); ++gg) {
141
142 // get integration weight and Jacobian of integration point (area of face)
143 double val = getGaussPts()(2, gg);
144 val *= cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
145 if (fe_type == MBTRI)
146 val /= 2;
147 for (int dd = 0; dd != 3; ++dd) {
148 coords[dd] = getCoordsAtGaussPts()(gg, dd);
149 normal[dd] = getNormalsAtGaussPts()(gg, dd);
150 }
151
152 if (analyticalForceOp) {
153 CHKERR analyticalForceOp->getForce(ent, coords, normal, force);
154 for (int rr = 0; rr != 3; ++rr) {
155 cblas_daxpy(nb_row_dofs, val * force[rr], &data.getN()(gg, 0), 1,
156 &nF[rr], rank);
157 }
158 } else {
159 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "No force to apply");
160 }
161 }
162
163 // Scale force using user defined scaling operator
165
166 auto get_f = [&]() {
167 if (F == PETSC_NULL)
168 return getKSPf();
169 return F;
170 };
171
172 // Assemble force into vector
173 CHKERR VecSetValues(get_f(), data, &*nF.data().begin(), ADD_VALUES);
174
176}
177
179 const std::string field_name, Vec _F, bCPressure &data,
180 boost::ptr_vector<MethodForForceScaling> &methods_op, bool ho_geometry)
183 F(_F), dAta(data), methodsOp(methods_op), hoGeometry(ho_geometry) {}
184
186 int side, EntityType type, EntitiesFieldData::EntData &data) {
187
189
190 if (data.getIndices().size() == 0)
192
193 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
194 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
195 if (dAta.tRis.find(fe_ent) == dAta.tRis.end())
197
198 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
199 int nb_row_dofs = data.getIndices().size() / rank;
200
201 Nf.resize(data.getIndices().size(), false);
202 Nf.clear();
203
204 for (unsigned int gg = 0; gg != data.getN().size1(); ++gg) {
205
206 double val = getGaussPts()(2, gg);
207 if (fe_type == MBTRI)
208 val /= 2;
209
210 for (int rr = 0; rr != rank; ++rr) {
211
212 double force;
213 if (hoGeometry || fe_type == MBQUAD)
214 force = dAta.data.data.value1 * getNormalsAtGaussPts()(gg, rr);
215 else
216 force = dAta.data.data.value1 * getNormal()[rr];
217
218 cblas_daxpy(nb_row_dofs, val * force, &data.getN()(gg, 0), 1, &Nf[rr],
219 rank);
220 }
221 }
222
224
225 auto get_f = [&]() {
226 if (F == PETSC_NULL)
227 return getKSPf();
228 return F;
229 };
230
231 CHKERR VecSetValues(get_f(), data, &*Nf.data().begin(), ADD_VALUES);
232
234}
235
237 int side, EntityType type, EntitiesFieldData::EntData &data) {
239
240 if (data.getFieldData().size() == 0)
241 PetscFunctionReturn(0);
242
243 ngp = data.getN().size1();
244
245 unsigned int nb_dofs = data.getFieldData().size() / 3;
247 if (type == MBVERTEX) {
248 dataAtIntegrationPts->tangent1.resize(3, ngp, false);
249 dataAtIntegrationPts->tangent1.clear();
250
251 dataAtIntegrationPts->tangent2.resize(3, ngp, false);
252 dataAtIntegrationPts->tangent2.clear();
253 }
254
255 auto t_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
256 auto t_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
257
258 for (unsigned int gg = 0; gg != ngp; ++gg) {
259 auto t_N = data.getFTensor1DiffN<2>(gg, 0);
260 auto t_dof = data.getFTensor1FieldData<3>();
261
262 for (unsigned int dd = 0; dd != nb_dofs; ++dd) {
263 t_1(i) += t_dof(i) * t_N(0);
264 t_2(i) += t_dof(i) * t_N(1);
265 ++t_dof;
266 ++t_N;
267 }
268 ++t_1;
269 ++t_2;
270 }
271
273}
274
276 int row_side, int col_side, EntityType row_type, EntityType col_type,
278 EntitiesFieldData::EntData &col_data) {
280
281 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
282 dAta.tRis.end()) {
284 }
285
286 const int row_nb_dofs = row_data.getIndices().size();
287 if (!row_nb_dofs)
289 const int col_nb_dofs = col_data.getIndices().size();
290 if (!col_nb_dofs)
292 const int nb_gauss_pts = row_data.getN().size1();
293
294 int nb_base_fun_row = row_data.getFieldData().size() / 3;
295 int nb_base_fun_col = col_data.getFieldData().size() / 3;
296
297 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
298 NN.clear();
299
303
304 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
306 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
307 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
308 &m(r + 2, c + 2));
309 };
310
311 auto make_vec_der = [](auto t_N, auto t_1, auto t_2) {
316 t_n(i, j) = 0;
317 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
318 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
319 return t_n;
320 };
321
322 double lambda = 1;
323 if (auto arc_dof = dataAtIntegrationPts->arcLengthDof.lock()) {
324 lambda = arc_dof->getFieldData();
325 }
326
327 auto t_w = getFTensor0IntegrationWeight();
328 auto t_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
329 auto t_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
330 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
331
332 double val = 0.5 * t_w * lambda;
333
334 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
335
336 int bbc = 0;
337 for (; bbc != nb_base_fun_col; ++bbc) {
338 auto t_base = row_data.getFTensor0N(gg, 0);
339
340 int bbr = 0;
341 for (; bbr != nb_base_fun_row; ++bbr) {
342
343 auto t_d_n = make_vec_der(t_N, t_1, t_2);
344
345 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
346
347 // TODO: handle hoGeometry
348 t_assemble(i, k) += val * dAta.data.data.value1 * t_base * t_d_n(i, k);
349
350 ++t_base;
351 }
352 ++t_N;
353 }
354 ++t_w;
355 ++t_1;
356 ++t_2;
357 }
358
359 // get pointer to first global index on row
360 const int *row_indices = &*row_data.getIndices().data().begin();
361 // get pointer to first global index on column
362 const int *col_indices = &*col_data.getIndices().data().begin();
363
364 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
365 : getFEMethod()->snes_B;
366
367 CHKERR MatSetValues(B, row_nb_dofs, row_indices, col_nb_dofs, col_indices,
368 &*NN.data().begin(), ADD_VALUES);
369
371}
372
374 int side, EntityType type, EntitiesFieldData::EntData &row_data) {
375
380
381 // get number of integration points
382 const int nb_integration_pts = getGaussPts().size2();
383 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hMat);
384 auto t_H = getFTensor2FromMat<3, 3>(dataAtPts->HMat);
385
386 dataAtPts->detHVec.resize(nb_integration_pts, false);
387 dataAtPts->invHMat.resize(9, nb_integration_pts, false);
388 dataAtPts->FMat.resize(9, nb_integration_pts, false);
389
390 auto t_detH = getFTensor0FromVec(dataAtPts->detHVec);
391 auto t_invH = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
392 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
393
394 for (int gg = 0; gg != nb_integration_pts; ++gg) {
395 CHKERR determinantTensor3by3(t_H, t_detH);
396 CHKERR invertTensor3by3(t_H, t_detH, t_invH);
397 t_F(i, j) = t_h(i, k) * t_invH(k, j);
398
399 ++t_h;
400 ++t_H;
401 ++t_detH;
402 ++t_invH;
403 ++t_F;
404 }
405
407}
408
410 int side, EntityType type, EntitiesFieldData::EntData &row_data) {
411
413
414 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
415 dAta.tRis.end()) {
417 }
418
419 // get number of dofs on row
420 nbRows = row_data.getIndices().size();
421 // if no dofs on row, exit that work, nothing to do here
422 if (!nbRows)
424
425 nF.resize(nbRows, false);
426 nF.clear();
427
428 // get number of integration points
429 nbIntegrationPts = getGaussPts().size2();
430
431 // integrate local matrix for entity block
432 CHKERR iNtegrate(row_data);
433
434 // assemble local matrix
435 CHKERR aSsemble(row_data);
436
438}
439
441 EntData &data) {
443
444 auto get_tensor1 = [](VectorDouble &v, const int r) {
446 &v(r + 0), &v(r + 1), &v(r + 2));
447 };
448
451
452 CHKERR loopSideVolumes(sideFeName, *sideFe);
453
454 // get integration weights
455 auto t_w = getFTensor0IntegrationWeight();
456
457 auto t_normal = getFTensor1NormalsAtGaussPts();
458
459 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
460
461 // iterate over integration points
462 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
463
464 FTensor::Tensor0<double *> t_base(&data.getN()(gg, 0));
465
466 double a = -0.5 * t_w * dAta.data.data.value1;
467 auto t_nf = get_tensor1(nF, 0);
468
469 int rr = 0;
470 for (; rr != nbRows / 3; ++rr) {
471
472 t_nf(i) += a * t_base * t_F(j, i) * t_normal(j);
473
474 ++t_base;
475 ++t_nf;
476 }
477
478 ++t_w;
479 ++t_F;
480 ++t_normal;
481 }
482
484}
485
487 EntData &row_data) {
489
490 // get pointer to first global index on row
491 const int *row_indices = &*row_data.getIndices().data().begin();
492
493 auto &data = *dataAtPts;
494 if (!data.forcesOnlyOnEntitiesRow.empty()) {
495 rowIndices.resize(nbRows, false);
496 noalias(rowIndices) = row_data.getIndices();
497 row_indices = &rowIndices[0];
498 VectorDofs &dofs = row_data.getFieldDofs();
499 VectorDofs::iterator dit = dofs.begin();
500 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
501 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
502 data.forcesOnlyOnEntitiesRow.end()) {
503 rowIndices[ii] = -1;
504 }
505 }
506 }
507
508 auto get_f = [&]() {
509 if (F == PETSC_NULL)
510 return getKSPf();
511 return F;
512 };
513
514 auto vec_assemble = [&](Vec my_f) {
516 CHKERR VecSetOption(my_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
517 CHKERR VecSetValues(my_f, nbRows, row_indices, &*nF.data().begin(),
518 ADD_VALUES);
520 };
521
522 // assemble local matrix
523 CHKERR vec_assemble(get_f());
524
526}
527
529 int row_side, int col_side, EntityType row_type, EntityType col_type,
531 EntitiesFieldData::EntData &col_data) {
532
534
535 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
536 dAta.tRis.end()) {
538 }
539
540 row_nb_dofs = row_data.getIndices().size();
541 if (!row_nb_dofs)
543 col_nb_dofs = col_data.getIndices().size();
544 if (!col_nb_dofs)
546 nb_gauss_pts = row_data.getN().size1();
547
548 nb_base_fun_row = row_data.getFieldData().size() / 3;
549 nb_base_fun_col = col_data.getFieldData().size() / 3;
550
551 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
552 NN.clear();
553
554 diagonal_block = (row_type == col_type) && (row_side == col_side);
555
556 if (col_type == MBVERTEX) {
557 dataAtPts->faceRowData = &row_data;
558 CHKERR loopSideVolumes(sideFeName, *sideFe);
559 }
560
561 // integrate local matrix for entity block
562 CHKERR iNtegrate(row_data, col_data);
563
564 // assemble local matrix
565 CHKERR aSsemble(row_data, col_data);
566
568}
569
572 EntData &row_data, EntData &col_data) {
573
575
579
580 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
582 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
583 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
584 &m(r + 2, c + 2));
585 };
586
587 auto make_vec_der = [](auto t_N, auto t_1, auto t_2) {
592 t_n(i, j) = 0;
593 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
594 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
595 return t_n;
596 };
597
598 dataAtPts->faceRowData = nullptr;
599 CHKERR loopSideVolumes(sideFeName, *sideFe);
600
601 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
602
603 double lambda = 1;
604 if (auto arc_dof = dataAtPts->arcLengthDof.lock()) {
605 lambda = arc_dof->getFieldData();
606 }
607
608 auto t_w = getFTensor0IntegrationWeight();
609 auto t_1 = getFTensor1FromMat<3>(dataAtPts->tangent1);
610 auto t_2 = getFTensor1FromMat<3>(dataAtPts->tangent2);
611
612 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
613
614 double val = 0.5 * t_w * lambda;
615
616 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
617
618 int bbc = 0;
619 for (; bbc != nb_base_fun_col; ++bbc) {
620
621 auto t_base = row_data.getFTensor0N(gg, 0);
622
623 int bbr = 0;
624 for (; bbr != nb_base_fun_row; ++bbr) {
625
626 auto t_d_n = make_vec_der(t_N, t_1, t_2);
627
628 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
629 // TODO: handle hoGeometry
630
631 t_assemble(i, k) -=
632 val * dAta.data.data.value1 * t_base * t_F(j, i) * t_d_n(j, k);
633
634 ++t_base;
635 }
636 ++t_N;
637 }
638 ++t_F;
639 ++t_w;
640 ++t_1;
641 ++t_2;
642 }
643
645}
646
648 EntData &row_data, EntData &col_data) {
649
651
652 // get pointer to first global index on row
653 const int *row_indices = &*row_data.getIndices().data().begin();
654 // get pointer to first global index on column
655 const int *col_indices = &*col_data.getIndices().data().begin();
656
657 auto &data = *dataAtPts;
658 if (!data.forcesOnlyOnEntitiesRow.empty()) {
659 rowIndices.resize(row_nb_dofs, false);
660 noalias(rowIndices) = row_data.getIndices();
661 row_indices = &rowIndices[0];
662 VectorDofs &dofs = row_data.getFieldDofs();
663 VectorDofs::iterator dit = dofs.begin();
664 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
665 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
666 data.forcesOnlyOnEntitiesRow.end()) {
667 rowIndices[ii] = -1;
668 }
669 }
670 }
671
672 if (!data.forcesOnlyOnEntitiesCol.empty()) {
673 colIndices.resize(col_nb_dofs, false);
674 noalias(colIndices) = col_data.getIndices();
675 col_indices = &colIndices[0];
676 VectorDofs &dofs = col_data.getFieldDofs();
677 VectorDofs::iterator dit = dofs.begin();
678 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
679 if (data.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
680 data.forcesOnlyOnEntitiesCol.end()) {
681 colIndices[ii] = -1;
682 }
683 }
684 }
685
686 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
687 : getFEMethod()->snes_B;
688 // assemble local matrix
689 CHKERR MatSetValues(B, row_nb_dofs, row_indices, col_nb_dofs, col_indices,
690 &*NN.data().begin(), ADD_VALUES);
691
693}
694
696 int row_side, int col_side, EntityType row_type, EntityType col_type,
698 EntitiesFieldData::EntData &col_data) {
699
701
702 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
703 dAta.tRis.end()) {
705 }
706
707 if (col_type != MBVERTEX)
709
710 dataAtPts->faceRowData = &row_data;
711 CHKERR loopSideVolumes(sideFeName, *sideFe);
712
714}
715
718 int row_side, int col_side, EntityType row_type, EntityType col_type,
720 EntitiesFieldData::EntData &col_data) {
721
723
724 if (dataAtPts->faceRowData == nullptr)
726
727 if (row_type != MBVERTEX)
729
730 row_nb_dofs = dataAtPts->faceRowData->getIndices().size();
731 if (!row_nb_dofs)
733 col_nb_dofs = col_data.getIndices().size();
734 if (!col_nb_dofs)
736
737 nb_gauss_pts = dataAtPts->faceRowData->getN().size1();
738
739 nb_base_fun_row = dataAtPts->faceRowData->getFieldData().size() / 3;
740 nb_base_fun_col = col_data.getFieldData().size() / 3;
741
742 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
743 NN.clear();
744
745 // integrate local matrix for entity block
746 CHKERR iNtegrate(*(dataAtPts->faceRowData), col_data);
747
748 // assemble local matrix
749 CHKERR aSsemble(*(dataAtPts->faceRowData), col_data);
750
752}
753
756 EntData &row_data, EntData &col_data) {
757
759
763
764 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
766 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
767 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
768 &m(r + 2, c + 2));
769 };
770
771 auto t_w = getFTensor0IntegrationWeight();
772
773 auto t_normal = getFTensor1NormalsAtGaussPts();
774
775 auto t_inv_H = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
776
777 double lambda = 1;
778 if (auto arc_dof = dataAtPts->arcLengthDof.lock()) {
779 lambda = arc_dof->getFieldData();
780 }
781
782 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
783
784 double a = -0.5 * t_w * dAta.data.data.value1 * lambda;
785
786 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
787
788 int bbc = 0;
789 for (; bbc != nb_base_fun_col; ++bbc) {
790
791 FTensor::Tensor0<double *> t_row_base(&row_data.getN()(gg, 0));
792
793 int bbr = 0;
794 for (; bbr != nb_base_fun_row; ++bbr) {
795
796 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
797 // TODO: handle hoGeometry
798
799 t_assemble(i, j) +=
800 a * t_row_base * t_inv_H(k, i) * t_col_diff_base(k) * t_normal(j);
801
802 ++t_row_base;
803 }
804 ++t_col_diff_base;
805 }
806 ++t_w;
807 ++t_normal;
808 ++t_inv_H;
809 }
810
812}
813
816 EntData &row_data, EntData &col_data) {
817
819
820 // get pointer to first global index on row
821 const int *row_indices = &*row_data.getIndices().data().begin();
822 // get pointer to first global index on column
823 const int *col_indices = &*col_data.getIndices().data().begin();
824
825 auto &data = *dataAtPts;
826 if (!data.forcesOnlyOnEntitiesRow.empty()) {
827 rowIndices.resize(row_nb_dofs, false);
828 noalias(rowIndices) = row_data.getIndices();
829 row_indices = &rowIndices[0];
830 VectorDofs &dofs = row_data.getFieldDofs();
831 VectorDofs::iterator dit = dofs.begin();
832 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
833 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
834 data.forcesOnlyOnEntitiesRow.end()) {
835 rowIndices[ii] = -1;
836 }
837 }
838 }
839
840 if (!data.forcesOnlyOnEntitiesCol.empty()) {
841 colIndices.resize(col_nb_dofs, false);
842 noalias(colIndices) = col_data.getIndices();
843 col_indices = &colIndices[0];
844 VectorDofs &dofs = col_data.getFieldDofs();
845 VectorDofs::iterator dit = dofs.begin();
846 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
847 if (data.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
848 data.forcesOnlyOnEntitiesCol.end()) {
849 colIndices[ii] = -1;
850 }
851 }
852 }
853
854 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
855 : getFEMethod()->snes_B;
856 // assemble local matrix
857 CHKERR MatSetValues(B, row_nb_dofs, row_indices, col_nb_dofs, col_indices,
858 &*NN.data().begin(), ADD_VALUES);
859
861}
862
865 EntData &row_data, EntData &col_data) {
866
868
874
875 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
877 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
878 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
879 &m(r + 2, c + 2));
880 };
881
882 auto t_w = getFTensor0IntegrationWeight();
883
884 auto t_normal = getFTensor1NormalsAtGaussPts();
885
886 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hMat);
887 auto t_inv_H = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
888
890
891 double lambda = 1;
892 if (auto arc_dof = dataAtPts->arcLengthDof.lock()) {
893 lambda = arc_dof->getFieldData();
894 }
895
896 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
897
898 double a = -0.5 * t_w * dAta.data.data.value1 * lambda;
899
900 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
901
902 int bbc = 0;
903 for (; bbc != nb_base_fun_col; ++bbc) {
904
905 FTensor::Tensor0<double *> t_row_base(&row_data.getN()(gg, 0));
906
907 int bbr = 0;
908 for (; bbr != nb_base_fun_row; ++bbr) {
909
910 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
911
912 // TODO: handle hoGeometry
913
914 t_assemble(i, j) -= a * t_row_base * t_inv_H(l, j) *
915 t_col_diff_base(m) * t_inv_H(m, i) * t_h(k, l) *
916 t_normal(k);
917
918 ++t_row_base;
919 }
920 ++t_col_diff_base;
921 }
922 ++t_w;
923 ++t_h;
924 ++t_inv_H;
925 ++t_normal;
926 }
927
929}
930
932 const std::string field_name, Vec _F, bCPressure &data,
933 boost::ptr_vector<MethodForForceScaling> &methods_op, bool ho_geometry)
936 F(_F), dAta(data), methodsOp(methods_op), hoGeometry(ho_geometry) {}
937
939 int side, EntityType type, EntitiesFieldData::EntData &data) {
941
942 if (data.getIndices().size() == 0)
944 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
945 dAta.tRis.end()) {
947 }
948
949 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
950 int nb_row_dofs = data.getIndices().size() / rank;
951
952 Nf.resize(data.getIndices().size(), false);
953 Nf.clear();
954
955 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
956
957 for (unsigned int gg = 0; gg != data.getN().size1(); ++gg) {
958
959 double val = getGaussPts()(2, gg);
960 double flux;
961 if (hoGeometry || fe_type == MBQUAD) {
962 double area = cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
963 if (fe_type == MBTRI)
964 area /= 2;
965 flux = dAta.data.data.value1 * area;
966 } else {
967 flux = dAta.data.data.value1 * getArea();
968 }
969
970 cblas_daxpy(nb_row_dofs, val * flux, &data.getN()(gg, 0), 1,
971 &*Nf.data().begin(), 1);
972 }
973
975
976 auto get_f = [&]() {
977 if (F == PETSC_NULL)
978 return getKSPf();
979 return F;
980 };
981
982 CHKERR VecSetValues(get_f(), data, &*Nf.data().begin(), ADD_VALUES);
983
985}
986
988 Vec F, int ms_id,
989 bool ho_geometry,
990 bool block_set) {
991 const CubitMeshSets *cubit_meshset_ptr;
992 MeshsetsManager *mmanager_ptr;
994 CHKERR mField.getInterface(mmanager_ptr);
995 if (block_set) {
996 // Add data from block set.
997 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, BLOCKSET,
998 &cubit_meshset_ptr);
999 std::vector<double> mydata;
1000 CHKERR cubit_meshset_ptr->getAttributes(mydata);
1001 VectorDouble force(mydata.size());
1002 for (unsigned int ii = 0; ii != mydata.size(); ++ii) {
1003 force[ii] = mydata[ii];
1004 }
1005 // Read forces from BLOCKSET Force (if exists)
1006 if (force.empty()) {
1007 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Force not given");
1008 }
1009 // Assign values from BLOCKSET FORCE to RHS vector. Info about native Cubit
1010 // BC data structure can be found in BCData.hpp
1011 const string name = "Force";
1012 strncpy(mapForce[ms_id].data.data.name, name.c_str(),
1013 name.size() > 5 ? 5 : name.size());
1014 double magnitude =
1015 std::sqrt(force[0] * force[0] + force[1] * force[1] + force[2] * force[2]);
1016 mapForce[ms_id].data.data.value1 = -magnitude; //< Force magnitude
1017 mapForce[ms_id].data.data.value2 = 0;
1018 mapForce[ms_id].data.data.value3 =
1019 force[0] / magnitude; //< X-component of force vector
1020 mapForce[ms_id].data.data.value4 =
1021 force[1] / magnitude; //< Y-component of force vector
1022 mapForce[ms_id].data.data.value5 =
1023 force[2] / magnitude; //< Z-component of force vector
1024 mapForce[ms_id].data.data.value6 = 0;
1025 mapForce[ms_id].data.data.value7 = 0;
1026 mapForce[ms_id].data.data.value8 = 0;
1027 mapForce[ms_id].data.data.zero[0] = 0;
1028 mapForce[ms_id].data.data.zero[1] = 0;
1029 mapForce[ms_id].data.data.zero[2] = 0;
1030 mapForce[ms_id].data.data.zero2 = 0;
1031
1032 CHKERR mField.get_moab().get_entities_by_dimension(
1033 cubit_meshset_ptr->meshset, 2, mapForce[ms_id].tRis, true);
1034 fe.getOpPtrVector().push_back(new OpNeumannForce(
1035 field_name, F, mapForce[ms_id], methodsOp, ho_geometry));
1036
1037 // SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,"Not implemented");
1038 } else {
1039 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, NODESET, &cubit_meshset_ptr);
1040 CHKERR cubit_meshset_ptr->getBcDataStructure(mapForce[ms_id].data);
1041 CHKERR mField.get_moab().get_entities_by_dimension(
1042 cubit_meshset_ptr->meshset, 2, mapForce[ms_id].tRis, true);
1043 fe.getOpPtrVector().push_back(new OpNeumannForce(
1044 field_name, F, mapForce[ms_id], methodsOp, ho_geometry));
1045 }
1047}
1048
1050 Vec F, int ms_id,
1051 bool ho_geometry,
1052 bool block_set) {
1053
1054 const CubitMeshSets *cubit_meshset_ptr;
1055 MeshsetsManager *mmanager_ptr;
1057 CHKERR mField.getInterface(mmanager_ptr);
1058 if (block_set) {
1059 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, BLOCKSET,
1060 &cubit_meshset_ptr);
1061 std::vector<double> mydata;
1062 CHKERR cubit_meshset_ptr->getAttributes(mydata);
1063 VectorDouble pressure(mydata.size());
1064 for (unsigned int ii = 0; ii != mydata.size(); ++ii) {
1065 pressure[ii] = mydata[ii];
1066 }
1067 if (pressure.empty()) {
1068 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Pressure not given");
1069 }
1070 const string name = "Pressure";
1071 strncpy(mapPressure[ms_id].data.data.name, name.c_str(),
1072 name.size() > 8 ? 8 : name.size());
1073 mapPressure[ms_id].data.data.flag1 = 0;
1074 mapPressure[ms_id].data.data.flag2 = 1;
1075 mapPressure[ms_id].data.data.value1 = pressure[0];
1076 mapPressure[ms_id].data.data.zero = 0;
1077 CHKERR mField.get_moab().get_entities_by_dimension(
1078 cubit_meshset_ptr->meshset, 2, mapPressure[ms_id].tRis, true);
1079 fe.getOpPtrVector().push_back(new OpNeumannPressure(
1080 field_name, F, mapPressure[ms_id], methodsOp, ho_geometry));
1081 } else {
1082 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, SIDESET, &cubit_meshset_ptr);
1083 CHKERR cubit_meshset_ptr->getBcDataStructure(mapPressure[ms_id].data);
1084 CHKERR mField.get_moab().get_entities_by_dimension(
1085 cubit_meshset_ptr->meshset, 2, mapPressure[ms_id].tRis, true);
1086 fe.getOpPtrVector().push_back(new OpNeumannPressure(
1087 field_name, F, mapPressure[ms_id], methodsOp, ho_geometry));
1088 }
1090}
1091
1093 const std::string x_field, const std::string X_field,
1094 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1095 std::string side_fe_name, Vec F, Mat aij, int ms_id, bool ho_geometry,
1096 bool block_set) {
1097
1098 const CubitMeshSets *cubit_meshset_ptr;
1099 MeshsetsManager *mmanager_ptr;
1101 CHKERR mField.getInterface(mmanager_ptr);
1102 if (block_set) {
1103 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, BLOCKSET,
1104 &cubit_meshset_ptr);
1105 std::vector<double> mydata;
1106 CHKERR cubit_meshset_ptr->getAttributes(mydata);
1107 VectorDouble pressure(mydata.size());
1108 for (unsigned int ii = 0; ii != mydata.size(); ++ii) {
1109 pressure[ii] = mydata[ii];
1110 }
1111 if (pressure.empty()) {
1112 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Pressure not given");
1113 }
1114 const string name = "Pressure";
1115 strncpy(mapPressure[ms_id].data.data.name, name.c_str(),
1116 name.size() > 8 ? 8 : name.size());
1117 mapPressure[ms_id].data.data.flag1 = 0;
1118 mapPressure[ms_id].data.data.flag2 = 1;
1119 mapPressure[ms_id].data.data.value1 = pressure[0];
1120 mapPressure[ms_id].data.data.zero = 0;
1121 CHKERR mField.get_moab().get_entities_by_dimension(
1122 cubit_meshset_ptr->meshset, 2, mapPressure[ms_id].tRis, true);
1123
1124 } else {
1125 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, SIDESET, &cubit_meshset_ptr);
1126 CHKERR cubit_meshset_ptr->getBcDataStructure(mapPressure[ms_id].data);
1127 CHKERR mField.get_moab().get_entities_by_dimension(
1128 cubit_meshset_ptr->meshset, 2, mapPressure[ms_id].tRis, true);
1129 }
1130
1131 /* LEFT-HAND SIDE (SPATIAL) */
1132
1133 feLhs.getOpPtrVector().push_back(new OpGetTangent(X_field, data_at_pts));
1134
1136 x_field, X_field, data_at_pts, aij, mapPressure[ms_id], ho_geometry));
1137
1138 /* RIGHT-HAND SIDE (MATERIAL) */
1139
1140 // Side volume element computes the deformation gradient F=hH^-1
1141 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideRhs =
1142 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1143
1144 feMatSideRhs->getOpPtrVector().push_back(
1146 data_at_pts->getHMatPtr()));
1147 feMatSideRhs->getOpPtrVector().push_back(
1149 data_at_pts->getSmallhMatPtr()));
1150
1151 feMatSideRhs->getOpPtrVector().push_back(
1152 new OpCalculateDeformation(X_field, data_at_pts, ho_geometry));
1153
1155 X_field, data_at_pts, feMatSideRhs, side_fe_name, F, mapPressure[ms_id],
1156 ho_geometry));
1157
1158 /* LEFT-HAND SIDE (MATERIAL) */
1159
1160 // Side volume element computes linearisation with spatial coordinates
1161 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dx =
1162 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1163 // Side volume element computes linearisation with material coordinates
1164 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dX =
1165 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1166
1167 feMatSideLhs_dx->getOpPtrVector().push_back(
1169 data_at_pts->getHMatPtr()));
1170 feMatSideLhs_dx->getOpPtrVector().push_back(
1172 data_at_pts->getSmallhMatPtr()));
1173
1174 feMatSideLhs_dx->getOpPtrVector().push_back(
1175 new OpCalculateDeformation(X_field, data_at_pts, ho_geometry));
1176 feMatSideLhs_dx->getOpPtrVector().push_back(
1178 X_field, x_field, data_at_pts, aij, mapPressure[ms_id], ho_geometry));
1179
1181 X_field, x_field, data_at_pts, feMatSideLhs_dx, side_fe_name, aij,
1182 mapPressure[ms_id], ho_geometry));
1183
1184 feMatSideLhs_dX->getOpPtrVector().push_back(
1186 data_at_pts->getHMatPtr()));
1187 feMatSideLhs_dX->getOpPtrVector().push_back(
1189 data_at_pts->getSmallhMatPtr()));
1190 feMatSideLhs_dX->getOpPtrVector().push_back(
1191 new OpCalculateDeformation(X_field, data_at_pts, ho_geometry));
1192 feMatSideLhs_dX->getOpPtrVector().push_back(
1194 X_field, X_field, data_at_pts, aij, mapPressure[ms_id], ho_geometry));
1195
1196 feMatLhs.getOpPtrVector().push_back(new OpGetTangent(X_field, data_at_pts));
1198 X_field, X_field, data_at_pts, feMatSideLhs_dX, side_fe_name, aij,
1199 mapPressure[ms_id], ho_geometry));
1200
1202}
1203
1206 int ms_id, bool ho_geometry) {
1207
1208 const CubitMeshSets *cubit_meshset_ptr;
1209 MeshsetsManager *mmanager_ptr;
1211 CHKERR mField.getInterface(mmanager_ptr);
1212 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, BLOCKSET, &cubit_meshset_ptr);
1213 std::vector<double> mydata;
1214 CHKERR cubit_meshset_ptr->getAttributes(mydata);
1215 if (mydata.size() < 4)
1216 SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1217 "Should be four block attributes but is %d", mydata.size());
1218 VectorDouble3 pressure_coeffs(3);
1219 for (unsigned int ii = 0; ii != 3; ++ii) {
1220 pressure_coeffs[ii] = mydata[ii];
1221 }
1222 const double pressure_shift = mydata[3];
1223
1224 Range tris;
1225 CHKERR mField.get_moab().get_entities_by_dimension(cubit_meshset_ptr->meshset,
1226 2, tris, true);
1227 boost::shared_ptr<MethodForAnalyticalForce> analytical_force_op(
1228 new LinearVaringPresssure(pressure_coeffs, pressure_shift));
1230 field_name, F, tris, methodsOp, analytical_force_op, ho_geometry));
1231
1233}
1234
1236 Vec F, int ms_id,
1237 bool ho_geometry,
1238 bool block_set) {
1239 return NeumannForcesSurface::addPressure(field_name, F, ms_id, ho_geometry,
1240 block_set);
1241}
1242
1244 Vec F, int ms_id,
1245 bool ho_geometry) {
1246 const CubitMeshSets *cubit_meshset_ptr;
1247 MeshsetsManager *mmanager_ptr;
1249 CHKERR mField.getInterface(mmanager_ptr);
1250 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, SIDESET, &cubit_meshset_ptr);
1251 CHKERR cubit_meshset_ptr->getBcDataStructure(mapPressure[ms_id].data);
1252 CHKERR mField.get_moab().get_entities_by_dimension(
1253 cubit_meshset_ptr->meshset, 2, mapPressure[ms_id].tRis, true);
1254 fe.getOpPtrVector().push_back(new OpNeumannFlux(
1255 field_name, F, mapPressure[ms_id], methodsOp, ho_geometry));
1257}
1258
1260 MoFEM::Interface &m_field, const std::string field_name,
1261 const std::string mesh_nodals_positions, Range *intersect_ptr) {
1263
1264 // Define boundary element that operates on rows, columns and data of a
1265 // given field
1266 CHKERR m_field.add_finite_element("FORCE_FE", MF_ZERO);
1270 if (m_field.check_field(mesh_nodals_positions)) {
1272 mesh_nodals_positions);
1273 }
1274 // Add entities to that element, here we add all triangles with FORCESET
1275 // from cubit
1277 it)) {
1278 Range tris;
1279 CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 2, tris,
1280 true);
1281 if (intersect_ptr)
1282 tris = intersect(tris, *intersect_ptr);
1283 CHKERR m_field.add_ents_to_finite_element_by_dim(tris, 2, "FORCE_FE");
1284 }
1285
1286 CHKERR m_field.add_finite_element("PRESSURE_FE", MF_ZERO);
1289 CHKERR m_field.modify_finite_element_add_field_data("PRESSURE_FE",
1290 field_name);
1291 if (m_field.check_field(mesh_nodals_positions)) {
1292 CHKERR m_field.modify_finite_element_add_field_data("PRESSURE_FE",
1293 mesh_nodals_positions);
1294 }
1295
1297 SIDESET | PRESSURESET, it)) {
1298 Range tris;
1299 CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 2, tris,
1300 true);
1301 if (intersect_ptr)
1302 tris = intersect(tris, *intersect_ptr);
1303 CHKERR m_field.add_ents_to_finite_element_by_dim(tris, 2, "PRESSURE_FE");
1304 }
1305
1306 // Reading forces from BLOCKSET
1307
1308 const string block_set_force_name("FORCE");
1309 // search for block named FORCE and add its attributes to FORCE_FE element
1311 if (it->getName().compare(0, block_set_force_name.length(),
1312 block_set_force_name) == 0) {
1313 Range tris;
1314 CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 2, tris,
1315 true);
1316 if (intersect_ptr)
1317 tris = intersect(tris, *intersect_ptr);
1318 CHKERR m_field.add_ents_to_finite_element_by_dim(tris, 2, "FORCE_FE");
1319 }
1320 }
1321 // search for block named PRESSURE and add its attributes to PRESSURE_FE
1322 // element
1323 const string block_set_pressure_name("PRESSURE");
1325 if (it->getName().compare(0, block_set_pressure_name.length(),
1326 block_set_pressure_name) == 0) {
1327 Range tris;
1328 CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 2, tris,
1329 true);
1330 if (intersect_ptr)
1331 tris = intersect(tris, *intersect_ptr);
1332 CHKERR m_field.add_ents_to_finite_element_by_dim(tris, 2, "PRESSURE_FE");
1333 }
1334 }
1335
1336 // search for block named LINEAR_PRESSURE and add its attributes to
1337 // PRESSURE_FE element
1338 const string block_set_linear_pressure_name("LINEAR_PRESSURE");
1340 if (it->getName().compare(0, block_set_linear_pressure_name.length(),
1341 block_set_linear_pressure_name) == 0) {
1342 Range tris;
1343 CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 2, tris,
1344 true);
1345 if (intersect_ptr)
1346 tris = intersect(tris, *intersect_ptr);
1347 CHKERR m_field.add_ents_to_finite_element_by_dim(tris, 2, "PRESSURE_FE");
1348 }
1349 }
1350
1352}
1353
1355 MoFEM::Interface &m_field,
1356 boost::ptr_map<std::string, NeumannForcesSurface> &neumann_forces, Vec F,
1357 const std::string field_name, const std::string mesh_nodals_positions) {
1359 bool ho_geometry = m_field.check_field(mesh_nodals_positions);
1360 // Add forces
1361 {
1362 std::string fe_name = "FORCE_FE";
1363 neumann_forces.insert(fe_name, new NeumannForcesSurface(m_field));
1364 auto &nf = neumann_forces.at(fe_name);
1365 auto &fe = nf.getLoopFe();
1366
1367 if (ho_geometry)
1368 CHKERR addHOOpsFace3D(mesh_nodals_positions, fe, false, false);
1369
1371 it)) {
1372 CHKERR nf.addForce(field_name, F, it->getMeshsetId(), ho_geometry, false);
1373 }
1374 // Reading forces from BLOCKSET
1375 const string block_set_force_name("FORCE");
1377 if (it->getName().compare(0, block_set_force_name.length(),
1378 block_set_force_name) == 0) {
1379 CHKERR nf.addForce(field_name, F, it->getMeshsetId(), ho_geometry,
1380 true);
1381 }
1382 }
1383 }
1384
1385 // Add pressures
1386 {
1387 std::string fe_name = "PRESSURE_FE";
1388 neumann_forces.insert(fe_name, new NeumannForcesSurface(m_field));
1389
1390 auto &nf = neumann_forces.at(fe_name);
1391 auto &fe = nf.getLoopFe();
1392
1393 if (ho_geometry)
1394 CHKERR addHOOpsFace3D(mesh_nodals_positions, fe, false, false);
1395
1397 m_field, SIDESET | PRESSURESET, it)) {
1398 CHKERR nf.addPressure(field_name, F, it->getMeshsetId(), ho_geometry,
1399 false);
1400 }
1401
1402 // Reading pressures from BLOCKSET
1403 const string block_set_pressure_name("PRESSURE");
1405 if (it->getName().compare(0, block_set_pressure_name.length(),
1406 block_set_pressure_name) == 0) {
1407 CHKERR nf.addPressure(field_name, F, it->getMeshsetId(), ho_geometry,
1408 true);
1409 }
1410 }
1411
1412 // Reading pressures from BLOCKSET
1413 const string block_set_linear_pressure_name("LINEAR_PRESSURE");
1415 if (it->getName().compare(0, block_set_linear_pressure_name.length(),
1416 block_set_linear_pressure_name) == 0) {
1417 CHKERR nf.addLinearPressure(field_name, F, it->getMeshsetId(),
1418 ho_geometry);
1419 }
1420 }
1421 }
1422
1424}
1425
1427 MoFEM::Interface &m_field, const std::string field_name,
1428 const std::string mesh_nodals_positions) {
1430
1431 CHKERR m_field.add_finite_element("FLUX_FE", MF_ZERO);
1435 if (m_field.check_field(mesh_nodals_positions)) {
1437 mesh_nodals_positions);
1438 }
1439
1441 SIDESET | PRESSURESET, it)) {
1442 Range tris;
1443 CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 2, tris,
1444 true);
1445 CHKERR m_field.add_ents_to_finite_element_by_dim(tris, 2, "FLUX_FE");
1446 }
1447
1449}
1450
1452 MoFEM::Interface &m_field,
1453 boost::ptr_map<std::string, NeumannForcesSurface> &neumann_forces, Vec F,
1454 const std::string field_name, const std::string mesh_nodals_positions) {
1456 bool ho_geometry = m_field.check_field(mesh_nodals_positions);
1457
1458 string fe_name;
1459 fe_name = "FLUX_FE";
1460 neumann_forces.insert(fe_name, new NeumannForcesSurface(m_field));
1461
1462 auto &nf = neumann_forces.at(fe_name);
1463 auto &fe = nf.getLoopFe();
1464
1465 if (ho_geometry)
1466 CHKERR addHOOpsFace3D(mesh_nodals_positions, fe, false, false);
1467
1469 SIDESET | PRESSURESET, it)) {
1470 CHKERR nf.addFlux(field_name, F, it->getMeshsetId(), ho_geometry);
1471 }
1473}
static Index< 'p', 3 > p
constexpr double a
@ MF_ZERO
Definition: definitions.h:111
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ PRESSURESET
Definition: definitions.h:165
@ FORCESET
Definition: definitions.h:164
@ NODESET
Definition: definitions.h:159
@ SIDESET
Definition: definitions.h:160
@ BLOCKSET
Definition: definitions.h:161
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_INVALID_DATA
Definition: definitions.h:49
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
constexpr double lambda
auto get_f
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
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
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
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
virtual bool check_field(const std::string &name) const =0
check if field is in database
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
#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.
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:77
const FTensor::Tensor2< T, Dim, Dim > Vec
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
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:103
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
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:1218
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1207
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
const double r
rate factor
double flux
impulse magnitude
constexpr auto field_name
#define _F(n)
Definition: quad.c:25
FTensor::Index< 'm', 3 > m
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")
static MoFEMErrorCode addNeumannFluxBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
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.
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.
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
virtual moab::Interface & get_moab()=0
this struct keeps basic methods for moab meshset about material and boundary conditions
MoFEMErrorCode getAttributes(std::vector< double > &attributes) const
get Cubit block attributes
MoFEMErrorCode getBcDataStructure(CUBIT_BC_DATA_TYPE &data) const
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Interface for managing meshsets containing materials and boundary conditions.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode getForce(const EntityHandle ent, const VectorDouble3 &coords, const VectorDouble3 &normal, VectorDouble3 &force)
MyTriangleFE(MoFEM::Interface &m_field)
Operator for computing deformation gradients in side volumes.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for computing tangent vectors.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpNeumannFlux(const std::string field_name, Vec _F, bCPressure &data, boost::ptr_vector< MethodForForceScaling > &methods_op, bool ho_geometry)
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)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpNeumannForce(const std::string field_name, Vec _F, bCForce &data, boost::ptr_vector< MethodForForceScaling > &methods_op, bool ho_geometry=false)
RHS-operator for pressure element (spatial configuration)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate pressure.
OpNeumannPressure(const std::string field_name, Vec _F, bCPressure &data, boost::ptr_vector< MethodForForceScaling > &methods_op, bool ho_geometry=false)
LHS-operator for pressure element (spatial configuration)
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.
LHS-operator for the pressure element (material configuration)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Compute part of the left-hand side.
LHS-operator for the pressure element (material configuration)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
RHS-operator for the pressure element (material configuration)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &row_data)
Integrate pressure in the material configuration.
LHS-operator (material configuration) on the side volume.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
LHS-operator (material configuration) on the side volume.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
std::map< int, bCPressure > mapPressure
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)
std::map< int, bCForce > mapForce
MoFEMErrorCode addFlux(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false)
Add flux element operator (integration on face)
boost::ptr_vector< MethodForForceScaling > methodsOp
DEPRECATED MoFEMErrorCode addPreassure(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false, bool block_set=false)
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.
NeumannForcesSurface(MoFEM::Interface &m_field)
MoFEMErrorCode addLinearPressure(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false)
Add operator to calculate pressure on element.
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.
boost::ptr_vector< MethodForAnalyticalForce > analyticalForceOp
MoFEM::Interface & mField