v0.13.1
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
15 const EntityHandle ent, const VectorDouble3 &coords,
16 const VectorDouble3 &normal, VectorDouble3 &force) {
18 const double p = inner_prod(coords, linearConstants) + pressureShift;
19 force = normal * p / norm_2(normal);
21}
22
24 : FaceElementForcesAndSourcesCore(m_field), addToRule(1) {}
25
27 const std::string field_name, Vec _F, bCForce &data,
28 boost::ptr_vector<MethodForForceScaling> &methods_op, bool ho_geometry)
31 F(_F), dAta(data), methodsOp(methods_op), hoGeometry(ho_geometry) {}
32
37
38 if (data.getIndices().size() == 0)
40 EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
41 if (dAta.tRis.find(ent) == dAta.tRis.end())
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
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
95}
96
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)
104 F(f), tRis(tris), methodsOp(methods_op),
105 analyticalForceOp(analytical_force_op), hoGeometry(ho_geometry) {}
106
108 int side, EntityType type, EntitiesFieldData::EntData &data) {
110
111 if (data.getIndices().size() == 0)
113 const EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
114 if (tRis.find(ent) == tRis.end())
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
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
165}
166
168 const std::string field_name, Vec _F, bCPressure &data,
169 boost::ptr_vector<MethodForForceScaling> &methods_op, bool ho_geometry)
172 F(_F), dAta(data), methodsOp(methods_op), hoGeometry(ho_geometry) {}
173
175 int side, EntityType type, EntitiesFieldData::EntData &data) {
176
178
179 if (data.getIndices().size() == 0)
181
182 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
183 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
184 if (dAta.tRis.find(fe_ent) == dAta.tRis.end())
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
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
223}
224
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;
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
263}
264
266 int row_side, int col_side, EntityType row_type, EntityType col_type,
268 EntitiesFieldData::EntData &col_data) {
270
271 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
272 dAta.tRis.end()) {
274 }
275
276 const int row_nb_dofs = row_data.getIndices().size();
277 if (!row_nb_dofs)
279 const int col_nb_dofs = col_data.getIndices().size();
280 if (!col_nb_dofs)
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
293
294 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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) {
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
361}
362
364 int row_side, int col_side, EntityType row_type, EntityType col_type,
366 EntitiesFieldData::EntData &col_data) {
368
369 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
370 dAta.tRis.end()) {
372 }
373
374 const int row_nb_dofs = row_data.getIndices().size();
375 if (!row_nb_dofs)
377 const int col_nb_dofs = col_data.getIndices().size();
378 if (!col_nb_dofs)
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
391
392 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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) {
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
468}
469
471 int side, EntityType type, EntitiesFieldData::EntData &row_data) {
472
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
504}
505
507 int side, EntityType type, EntitiesFieldData::EntData &row_data) {
508
510
511 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
512 dAta.tRis.end()) {
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)
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
535}
536
538 EntData &data) {
540
541 auto get_tensor1 = [](VectorDouble &v, const int r) {
543 &v(r + 0), &v(r + 1), &v(r + 2));
544 };
545
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
581}
582
584 EntData &row_data) {
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) {
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);
616 };
617
618 // assemble local matrix
619 CHKERR vec_assemble(get_f());
620
622}
623
626 int side, EntityType type, EntitiesFieldData::EntData &row_data) {
627
629
630 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
631 dAta.tRis.end()) {
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)
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
654}
655
658 EntData &data) {
660
661 auto get_tensor1 = [](VectorDouble &v, const int r) {
663 &v(r + 0), &v(r + 1), &v(r + 2));
664 };
665
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
706}
707
710 EntData &row_data) {
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) {
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);
742 };
743
744 // assemble local matrix
745 CHKERR vec_assemble(get_f());
746
748}
749
751 int row_side, int col_side, EntityType row_type, EntityType col_type,
753 EntitiesFieldData::EntData &col_data) {
754
756
757 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
758 dAta.tRis.end()) {
760 }
761
762 row_nb_dofs = row_data.getIndices().size();
763 if (!row_nb_dofs)
765 col_nb_dofs = col_data.getIndices().size();
766 if (!col_nb_dofs)
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
790}
791
794 EntData &row_data, EntData &col_data) {
795
797
801
802 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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) {
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
867}
868
870 EntData &row_data, EntData &col_data) {
871
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
900}
901
904 int row_side, int col_side, EntityType row_type, EntityType col_type,
906 EntitiesFieldData::EntData &col_data) {
907
909
910 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
911 dAta.tRis.end()) {
913 }
914
915 row_nb_dofs = row_data.getIndices().size();
916 if (!row_nb_dofs)
918 col_nb_dofs = col_data.getIndices().size();
919 if (!col_nb_dofs)
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
943}
944
947 EntData &row_data, EntData &col_data) {
948
950
955
956 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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) {
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
1028}
1029
1031 EntData &row_data, EntData &col_data) {
1032
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
1061}
1062
1064 int row_side, int col_side, EntityType row_type, EntityType col_type,
1066 EntitiesFieldData::EntData &col_data) {
1067
1069
1070 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1071 dAta.tRis.end()) {
1073 }
1074
1075 if (col_type != MBVERTEX)
1077
1078 dataAtPts->faceRowData = &row_data;
1079 CHKERR loopSideVolumes(sideFeName, *sideFe);
1080
1082}
1083
1086 int row_side, int col_side, EntityType row_type, EntityType col_type,
1088 EntitiesFieldData::EntData &col_data) {
1089
1091
1092 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1093 dAta.tRis.end()) {
1095 }
1096
1097 if (col_type != MBVERTEX)
1099
1100 dataAtPts->faceRowData = &row_data;
1101 CHKERR loopSideVolumes(sideFeName, *sideFe);
1102
1104}
1105
1108 int row_side, int col_side, EntityType row_type, EntityType col_type,
1110 EntitiesFieldData::EntData &col_data) {
1111
1113
1114 if (dataAtPts->faceRowData == nullptr)
1116
1117 if (row_type != MBVERTEX)
1119
1120 row_nb_dofs = dataAtPts->faceRowData->getIndices().size();
1121 if (!row_nb_dofs)
1123 col_nb_dofs = col_data.getIndices().size();
1124 if (!col_nb_dofs)
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
1142}
1143
1146 int row_side, int col_side, EntityType row_type, EntityType col_type,
1148 EntitiesFieldData::EntData &col_data) {
1149
1151
1152 if (dataAtPts->faceRowData == nullptr)
1154
1155 if (row_type != MBVERTEX)
1157
1158 row_nb_dofs = dataAtPts->faceRowData->getIndices().size();
1159 if (!row_nb_dofs)
1161 col_nb_dofs = col_data.getIndices().size();
1162 if (!col_nb_dofs)
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
1180}
1181
1184 EntData &row_data, EntData &col_data) {
1185
1187
1191
1192 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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
1240}
1241
1244 iNtegrate(EntData &row_data, EntData &col_data) {
1245
1247
1251
1252 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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
1305}
1306
1309 EntData &row_data, EntData &col_data) {
1310
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
1338}
1339
1342 EntData &row_data, EntData &col_data) {
1343
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
1373}
1374
1377 EntData &row_data, EntData &col_data) {
1378
1380
1386
1387 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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
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
1441}
1442
1445 iNtegrate(EntData &row_data, EntData &col_data) {
1446
1448
1454
1455 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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
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
1512}
1513
1515 const std::string field_name, Vec _F, bCPressure &data,
1516 boost::ptr_vector<MethodForForceScaling> &methods_op, bool ho_geometry)
1519 F(_F), dAta(data), methodsOp(methods_op), hoGeometry(ho_geometry) {}
1520
1525
1526 if (data.getIndices().size() == 0)
1528 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1529 dAta.tRis.end()) {
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
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
1569}
1570
1572 Vec F, int ms_id,
1573 bool ho_geometry,
1574 bool block_set) {
1575 const CubitMeshSets *cubit_meshset_ptr;
1576 MeshsetsManager *mmanager_ptr;
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 }
1631}
1632
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;
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
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(
1703 data_at_pts->getHMatPtr()));
1704 feMatSideRhs->getOpPtrVector().push_back(
1706 x_field, data_at_pts->getSmallhMatPtr()));
1707
1708 feMatSideRhs->getOpPtrVector().push_back(
1709 new OpCalculateDeformation(X_field, data_at_pts, ho_geometry));
1710
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(
1726 data_at_pts->getHMatPtr()));
1727 feMatSideLhs_dx->getOpPtrVector().push_back(
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(
1735 X_field, x_field, data_at_pts, aij, mapForce[ms_id], ho_geometry));
1736
1737 feMatLhs.getOpPtrVector().push_back(
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(
1744 data_at_pts->getHMatPtr()));
1745 feMatSideLhs_dX->getOpPtrVector().push_back(
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(
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(
1757 X_field, X_field, data_at_pts, feMatSideLhs_dX, side_fe_name, aij,
1758 mapForce[ms_id], ho_geometry));
1759 }
1760
1762}
1763
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;
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 }
1805}
1806
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;
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
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(
1861 data_at_pts->getHMatPtr()));
1862 feMatSideRhs->getOpPtrVector().push_back(
1864 data_at_pts->getSmallhMatPtr()));
1865
1866 feMatSideRhs->getOpPtrVector().push_back(
1867 new OpCalculateDeformation(X_field, data_at_pts, ho_geometry));
1868
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(
1884 data_at_pts->getHMatPtr()));
1885 feMatSideLhs_dx->getOpPtrVector().push_back(
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(
1893 X_field, x_field, data_at_pts, aij, mapPressure[ms_id], ho_geometry));
1894
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(
1901 data_at_pts->getHMatPtr()));
1902 feMatSideLhs_dX->getOpPtrVector().push_back(
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(
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));
1913 X_field, X_field, data_at_pts, feMatSideLhs_dX, side_fe_name, aij,
1914 mapPressure[ms_id], ho_geometry));
1915
1917}
1918
1921 int ms_id, bool ho_geometry) {
1922
1923 const CubitMeshSets *cubit_meshset_ptr;
1924 MeshsetsManager *mmanager_ptr;
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));
1945 field_name, F, tris, methodsOp, analytical_force_op, ho_geometry));
1946
1948}
1949
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
1959 Vec F, int ms_id,
1960 bool ho_geometry) {
1961 const CubitMeshSets *cubit_meshset_ptr;
1962 MeshsetsManager *mmanager_ptr;
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));
1972}
1973
1975 MoFEM::Interface &m_field, const std::string field_name,
1976 const std::string mesh_nodals_positions, Range *intersect_ptr) {
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);
1985 if (m_field.check_field(mesh_nodals_positions)) {
1987 mesh_nodals_positions);
1988 }
1989 // Add entities to that element, here we add all triangles with FORCESET
1990 // from cubit
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);
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
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
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");
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");
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
2067}
2068
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) {
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
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");
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
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");
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");
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
2139}
2140
2142 MoFEM::Interface &m_field, const std::string field_name,
2143 const std::string mesh_nodals_positions) {
2145
2146 CHKERR m_field.add_finite_element("FLUX_FE", MF_ZERO);
2150 if (m_field.check_field(mesh_nodals_positions)) {
2152 mesh_nodals_positions);
2153 }
2154
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
2164}
2165
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) {
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
2184 SIDESET | PRESSURESET, it)) {
2185 CHKERR nf.addFlux(field_name, F, it->getMeshsetId(), ho_geometry);
2186 }
2188}
static Index< 'p', 3 > p
constexpr double a
T data[Tensor_Dim]
@ MF_ZERO
Definition: definitions.h:98
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ PRESSURESET
Definition: definitions.h:152
@ FORCESET
Definition: definitions.h:151
@ NODESET
Definition: definitions.h:146
@ SIDESET
Definition: definitions.h:147
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Index< 'm', SPACE_DIM > m
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
const double c
speed of light (cm/ns)
const 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
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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:1323
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1312
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
constexpr auto field_name
#define _F(n)
Definition: quad.c:25
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)
Integrate surface force (traction)
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)
LHS-operator for surface force 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 surface force 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 surface force 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 surface force element (material configuration)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &row_data)
Integrate surface force 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 aSsemble(EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Finite element and operators to apply force/pressures applied to surfaces.
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)
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)
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.
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