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
26MoFEMErrorCode NeumannForcesSurface::LinearVaringPresssure::getForce(
27 const EntityHandle ent, const VectorDouble3 &coords,
28 const VectorDouble3 &normal, VectorDouble3 &force) {
29 MoFEMFunctionBegin;
30 const double p = inner_prod(coords, linearConstants) + pressureShift;
31 force = normal * p / norm_2(normal);
32 MoFEMFunctionReturn(0);
33}
34
35NeumannForcesSurface::MyTriangleFE::MyTriangleFE(MoFEM::Interface &m_field)
36 : FaceElementForcesAndSourcesCore(m_field), addToRule(1) {}
37
38NeumannForcesSurface::OpNeumannForce::OpNeumannForce(
39 const std::string field_name, Vec _F, bCForce &data,
40 boost::ptr_vector<MethodForForceScaling> &methods_op, bool ho_geometry)
41 : FaceElementForcesAndSourcesCore::UserDataOperator(
42 field_name, UserDataOperator::OPROW),
43 F(_F), dAta(data), methodsOp(methods_op), hoGeometry(ho_geometry) {}
44
45MoFEMErrorCode NeumannForcesSurface::OpNeumannForce::doWork(
46 int side, EntityType type, EntitiesFieldData::EntData &data) {
47 MoFEMFunctionBegin;
48
49 if (data.getIndices().size() == 0)
50 MoFEMFunctionReturnHot(0);
51 EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
52 if (dAta.tRis.find(ent) == dAta.tRis.end())
53 MoFEMFunctionReturnHot(0);
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
94 CHKERR MethodForForceScaling::applyScale(getFEMethod(), methodsOp, Nf);
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
105 MoFEMFunctionReturn(0);
106}
107
108NeumannForcesSurface::OpNeumannForceAnalytical::OpNeumannForceAnalytical(
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)
113 : FaceElementForcesAndSourcesCore::UserDataOperator(
114 field_name, UserDataOperator::OPROW),
115 F(f), tRis(tris), methodsOp(methods_op),
116 analyticalForceOp(analytical_force_op), hoGeometry(ho_geometry) {}
117
118MoFEMErrorCode NeumannForcesSurface::OpNeumannForceAnalytical::doWork(
119 int side, EntityType type, EntitiesFieldData::EntData &data) {
120 MoFEMFunctionBegin;
121
122 if (data.getIndices().size() == 0)
123 MoFEMFunctionReturnHot(0);
124 const EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
125 if (tRis.find(ent) == tRis.end())
126 MoFEMFunctionReturnHot(0);
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);
137 VectorDouble3 normal(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
164 CHKERR MethodForForceScaling::applyScale(getFEMethod(), methodsOp, nF);
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
175 MoFEMFunctionReturn(0);
176}
177
178NeumannForcesSurface::OpNeumannPressure::OpNeumannPressure(
179 const std::string field_name, Vec _F, bCPressure &data,
180 boost::ptr_vector<MethodForForceScaling> &methods_op, bool ho_geometry)
181 : FaceElementForcesAndSourcesCore::UserDataOperator(
182 field_name, UserDataOperator::OPROW),
183 F(_F), dAta(data), methodsOp(methods_op), hoGeometry(ho_geometry) {}
184
185MoFEMErrorCode NeumannForcesSurface::OpNeumannPressure::doWork(
186 int side, EntityType type, EntitiesFieldData::EntData &data) {
187
188 MoFEMFunctionBegin;
189
190 if (data.getIndices().size() == 0)
191 MoFEMFunctionReturnHot(0);
192
193 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
194 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
195 if (dAta.tRis.find(fe_ent) == dAta.tRis.end())
196 MoFEMFunctionReturnHot(0);
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
223 CHKERR MethodForForceScaling::applyScale(getFEMethod(), methodsOp, Nf);
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
233 MoFEMFunctionReturn(0);
234}
235
236MoFEMErrorCode NeumannForcesSurface::OpGetTangent::doWork(
237 int side, EntityType type, EntitiesFieldData::EntData &data) {
238 MoFEMFunctionBegin;
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;
246 FTensor::Index<'i', 3> i;
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
272 MoFEMFunctionReturn(0);
273}
274
275MoFEMErrorCode NeumannForcesSurface::OpNeumannPressureLhs_dx_dX::doWork(
276 int row_side, int col_side, EntityType row_type, EntityType col_type,
277 EntitiesFieldData::EntData &row_data,
278 EntitiesFieldData::EntData &col_data) {
279 MoFEMFunctionBegin;
280
281 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
282 dAta.tRis.end()) {
283 MoFEMFunctionReturnHot(0);
284 }
285
286 const int row_nb_dofs = row_data.getIndices().size();
287 if (!row_nb_dofs)
288 MoFEMFunctionReturnHot(0);
289 const int col_nb_dofs = col_data.getIndices().size();
290 if (!col_nb_dofs)
291 MoFEMFunctionReturnHot(0);
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
300 FTensor::Index<'i', 3> i;
301 FTensor::Index<'j', 3> j;
302 FTensor::Index<'k', 3> k;
303
304 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
305 return FTensor::Tensor2<double *, 3, 3>(
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) {
312 FTensor::Index<'i', 3> i;
313 FTensor::Index<'j', 3> j;
314 FTensor::Index<'k', 3> k;
315 FTensor::Tensor2<double, 3, 3> t_n;
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
370 MoFEMFunctionReturn(0);
371}
372
373MoFEMErrorCode NeumannForcesSurface::OpCalculateDeformation::doWork(
374 int side, EntityType type, EntitiesFieldData::EntData &row_data) {
375
376 MoFEMFunctionBegin;
377 FTensor::Index<'i', 3> i;
378 FTensor::Index<'j', 3> j;
379 FTensor::Index<'k', 3> k;
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
406 MoFEMFunctionReturn(0);
407}
408
409MoFEMErrorCode NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX::doWork(
410 int side, EntityType type, EntitiesFieldData::EntData &row_data) {
411
412 MoFEMFunctionBegin;
413
414 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
415 dAta.tRis.end()) {
416 MoFEMFunctionReturnHot(0);
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)
423 MoFEMFunctionReturnHot(0);
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
437 MoFEMFunctionReturn(0);
438}
439
440MoFEMErrorCode NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX::iNtegrate(
441 EntData &data) {
442 MoFEMFunctionBegin;
443
444 auto get_tensor1 = [](VectorDouble &v, const int r) {
445 return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
446 &v(r + 0), &v(r + 1), &v(r + 2));
447 };
448
449 FTensor::Index<'i', 3> i;
450 FTensor::Index<'j', 3> j;
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
483 MoFEMFunctionReturn(0);
484}
485
486MoFEMErrorCode NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX::aSsemble(
487 EntData &row_data) {
488 MoFEMFunctionBegin;
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) {
515 MoFEMFunctionBegin;
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);
519 MoFEMFunctionReturn(0);
520 };
521
522 // assemble local matrix
523 CHKERR vec_assemble(get_f());
524
525 MoFEMFunctionReturn(0);
526}
527
528MoFEMErrorCode NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dX::doWork(
529 int row_side, int col_side, EntityType row_type, EntityType col_type,
530 EntitiesFieldData::EntData &row_data,
531 EntitiesFieldData::EntData &col_data) {
532
533 MoFEMFunctionBegin;
534
535 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
536 dAta.tRis.end()) {
537 MoFEMFunctionReturnHot(0);
538 }
539
540 row_nb_dofs = row_data.getIndices().size();
541 if (!row_nb_dofs)
542 MoFEMFunctionReturnHot(0);
543 col_nb_dofs = col_data.getIndices().size();
544 if (!col_nb_dofs)
545 MoFEMFunctionReturnHot(0);
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
567 MoFEMFunctionReturn(0);
568}
569
570MoFEMErrorCode
571NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dX::iNtegrate(
572 EntData &row_data, EntData &col_data) {
573
574 MoFEMFunctionBegin;
575
576 FTensor::Index<'i', 3> i;
577 FTensor::Index<'j', 3> j;
578 FTensor::Index<'k', 3> k;
579
580 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
581 return FTensor::Tensor2<double *, 3, 3>(
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) {
588 FTensor::Index<'i', 3> i;
589 FTensor::Index<'j', 3> j;
590 FTensor::Index<'k', 3> k;
591 FTensor::Tensor2<double, 3, 3> t_n;
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
644 MoFEMFunctionReturn(0);
645}
646
647MoFEMErrorCode NeumannForcesSurface::OpNeumannPressureMaterialLhs::aSsemble(
648 EntData &row_data, EntData &col_data) {
649
650 MoFEMFunctionBegin;
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
692 MoFEMFunctionReturn(0);
693}
694
695MoFEMErrorCode NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dx::doWork(
696 int row_side, int col_side, EntityType row_type, EntityType col_type,
697 EntitiesFieldData::EntData &row_data,
698 EntitiesFieldData::EntData &col_data) {
699
700 MoFEMFunctionBegin;
701
702 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
703 dAta.tRis.end()) {
704 MoFEMFunctionReturnHot(0);
705 }
706
707 if (col_type != MBVERTEX)
708 MoFEMFunctionReturnHot(0);
709
710 dataAtPts->faceRowData = &row_data;
711 CHKERR loopSideVolumes(sideFeName, *sideFe);
712
713 MoFEMFunctionReturn(0);
714}
715
716MoFEMErrorCode
717NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs::doWork(
718 int row_side, int col_side, EntityType row_type, EntityType col_type,
719 EntitiesFieldData::EntData &row_data,
720 EntitiesFieldData::EntData &col_data) {
721
722 MoFEMFunctionBegin;
723
724 if (dataAtPts->faceRowData == nullptr)
725 MoFEMFunctionReturnHot(0);
726
727 if (row_type != MBVERTEX)
728 MoFEMFunctionReturnHot(0);
729
730 row_nb_dofs = dataAtPts->faceRowData->getIndices().size();
731 if (!row_nb_dofs)
732 MoFEMFunctionReturnHot(0);
733 col_nb_dofs = col_data.getIndices().size();
734 if (!col_nb_dofs)
735 MoFEMFunctionReturnHot(0);
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
751 MoFEMFunctionReturn(0);
752}
753
754MoFEMErrorCode
755NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs_dX_dx::iNtegrate(
756 EntData &row_data, EntData &col_data) {
757
758 MoFEMFunctionBegin;
759
760 FTensor::Index<'i', 3> i;
761 FTensor::Index<'j', 3> j;
762 FTensor::Index<'k', 3> k;
763
764 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
765 return FTensor::Tensor2<double *, 3, 3>(
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
811 MoFEMFunctionReturn(0);
812}
813
814MoFEMErrorCode
815NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs::aSsemble(
816 EntData &row_data, EntData &col_data) {
817
818 MoFEMFunctionBegin;
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
860 MoFEMFunctionReturn(0);
861}
862
863MoFEMErrorCode
864NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs_dX_dX::iNtegrate(
865 EntData &row_data, EntData &col_data) {
866
867 MoFEMFunctionBegin;
868
869 FTensor::Index<'i', 3> i;
870 FTensor::Index<'j', 3> j;
871 FTensor::Index<'k', 3> k;
872 FTensor::Index<'l', 3> l;
873 FTensor::Index<'m', 3> m;
874
875 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
876 return FTensor::Tensor2<double *, 3, 3>(
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
889 FTensor::Tensor2<double, 3, 3> t_d;
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
928 MoFEMFunctionReturn(0);
929}
930
931NeumannForcesSurface::OpNeumannFlux::OpNeumannFlux(
932 const std::string field_name, Vec _F, bCPressure &data,
933 boost::ptr_vector<MethodForForceScaling> &methods_op, bool ho_geometry)
934 : FaceElementForcesAndSourcesCore::UserDataOperator(
935 field_name, UserDataOperator::OPROW),
936 F(_F), dAta(data), methodsOp(methods_op), hoGeometry(ho_geometry) {}
937
938MoFEMErrorCode NeumannForcesSurface::OpNeumannFlux::doWork(
939 int side, EntityType type, EntitiesFieldData::EntData &data) {
940 MoFEMFunctionBegin;
941
942 if (data.getIndices().size() == 0)
943 MoFEMFunctionReturnHot(0);
944 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
945 dAta.tRis.end()) {
946 MoFEMFunctionReturnHot(0);
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
974 CHKERR MethodForForceScaling::applyScale(getFEMethod(), methodsOp, Nf);
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
984 MoFEMFunctionReturn(0);
985}
986
987MoFEMErrorCode NeumannForcesSurface::addForce(const std::string field_name,
988 Vec F, int ms_id,
989 bool ho_geometry,
990 bool block_set) {
991 const CubitMeshSets *cubit_meshset_ptr;
992 MeshsetsManager *mmanager_ptr;
993 MoFEMFunctionBegin;
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 }
1046 MoFEMFunctionReturn(0);
1047}
1048
1049MoFEMErrorCode NeumannForcesSurface::addPressure(const std::string field_name,
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;
1056 MoFEMFunctionBegin;
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 }
1089 MoFEMFunctionReturn(0);
1090}
1091
1092MoFEMErrorCode NeumannForcesSurface::addPressureAle(
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;
1100 MoFEMFunctionBegin;
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
1135 feLhs.getOpPtrVector().push_back(new OpNeumannPressureLhs_dx_dX(
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(
1145 new OpCalculateVectorFieldGradient<3, 3>(X_field,
1146 data_at_pts->getHMatPtr()));
1147 feMatSideRhs->getOpPtrVector().push_back(
1148 new OpCalculateVectorFieldGradient<3, 3>(x_field,
1149 data_at_pts->getSmallhMatPtr()));
1150
1151 feMatSideRhs->getOpPtrVector().push_back(
1152 new OpCalculateDeformation(X_field, data_at_pts, ho_geometry));
1153
1154 feMatRhs.getOpPtrVector().push_back(new OpNeumannPressureMaterialRhs_dX(
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(
1168 new OpCalculateVectorFieldGradient<3, 3>(X_field,
1169 data_at_pts->getHMatPtr()));
1170 feMatSideLhs_dx->getOpPtrVector().push_back(
1171 new OpCalculateVectorFieldGradient<3, 3>(x_field,
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(
1177 new OpNeumannPressureMaterialVolOnSideLhs_dX_dx(
1178 X_field, x_field, data_at_pts, aij, mapPressure[ms_id], ho_geometry));
1179
1180 feMatLhs.getOpPtrVector().push_back(new OpNeumannPressureMaterialLhs_dX_dx(
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(
1185 new OpCalculateVectorFieldGradient<3, 3>(X_field,
1186 data_at_pts->getHMatPtr()));
1187 feMatSideLhs_dX->getOpPtrVector().push_back(
1188 new OpCalculateVectorFieldGradient<3, 3>(x_field,
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(
1193 new OpNeumannPressureMaterialVolOnSideLhs_dX_dX(
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));
1197 feMatLhs.getOpPtrVector().push_back(new OpNeumannPressureMaterialLhs_dX_dX(
1198 X_field, X_field, data_at_pts, feMatSideLhs_dX, side_fe_name, aij,
1199 mapPressure[ms_id], ho_geometry));
1200
1201 MoFEMFunctionReturn(0);
1202}
1203
1204MoFEMErrorCode
1205NeumannForcesSurface::addLinearPressure(const std::string field_name, Vec F,
1206 int ms_id, bool ho_geometry) {
1207
1208 const CubitMeshSets *cubit_meshset_ptr;
1209 MeshsetsManager *mmanager_ptr;
1210 MoFEMFunctionBegin;
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));
1229 fe.getOpPtrVector().push_back(new OpNeumannForceAnalytical(
1230 field_name, F, tris, methodsOp, analytical_force_op, ho_geometry));
1231
1232 MoFEMFunctionReturn(0);
1233}
1234
1235MoFEMErrorCode NeumannForcesSurface::addPreassure(const std::string field_name,
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
1243MoFEMErrorCode NeumannForcesSurface::addFlux(const std::string field_name,
1244 Vec F, int ms_id,
1245 bool ho_geometry) {
1246 const CubitMeshSets *cubit_meshset_ptr;
1247 MeshsetsManager *mmanager_ptr;
1248 MoFEMFunctionBegin;
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));
1256 MoFEMFunctionReturn(0);
1257}
1258
1259MoFEMErrorCode MetaNeumannForces::addNeumannBCElements(
1260 MoFEM::Interface &m_field, const std::string field_name,
1261 const std::string mesh_nodals_positions, Range *intersect_ptr) {
1262 MoFEMFunctionBegin;
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);
1267 CHKERR m_field.modify_finite_element_add_field_row("FORCE_FE", field_name);
1268 CHKERR m_field.modify_finite_element_add_field_col("FORCE_FE", field_name);
1269 CHKERR m_field.modify_finite_element_add_field_data("FORCE_FE", field_name);
1270 if (m_field.check_field(mesh_nodals_positions)) {
1271 CHKERR m_field.modify_finite_element_add_field_data("FORCE_FE",
1272 mesh_nodals_positions);
1273 }
1274 // Add entities to that element, here we add all triangles with FORCESET
1275 // from cubit
1276 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
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);
1287 CHKERR m_field.modify_finite_element_add_field_row("PRESSURE_FE", field_name);
1288 CHKERR m_field.modify_finite_element_add_field_col("PRESSURE_FE", field_name);
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
1296 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,
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
1310 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
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");
1324 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
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");
1339 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
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
1351 MoFEMFunctionReturn(0);
1352}
1353
1354MoFEMErrorCode MetaNeumannForces::setMomentumFluxOperators(
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) {
1358 MoFEMFunctionBegin;
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
1370 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
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");
1376 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
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
1396 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
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");
1404 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
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");
1414 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
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
1423 MoFEMFunctionReturn(0);
1424}
1425
1426MoFEMErrorCode MetaNeumannForces::addNeumannFluxBCElements(
1427 MoFEM::Interface &m_field, const std::string field_name,
1428 const std::string mesh_nodals_positions) {
1429 MoFEMFunctionBegin;
1430
1431 CHKERR m_field.add_finite_element("FLUX_FE", MF_ZERO);
1432 CHKERR m_field.modify_finite_element_add_field_row("FLUX_FE", field_name);
1433 CHKERR m_field.modify_finite_element_add_field_col("FLUX_FE", field_name);
1434 CHKERR m_field.modify_finite_element_add_field_data("FLUX_FE", field_name);
1435 if (m_field.check_field(mesh_nodals_positions)) {
1436 CHKERR m_field.modify_finite_element_add_field_data("FLUX_FE",
1437 mesh_nodals_positions);
1438 }
1439
1440 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,
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
1448 MoFEMFunctionReturn(0);
1449}
1450
1451MoFEMErrorCode MetaNeumannForces::setMassFluxOperators(
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) {
1455 MoFEMFunctionBegin;
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
1468 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,
1469 SIDESET | PRESSURESET, it)) {
1470 CHKERR nf.addFlux(field_name, F, it->getMeshsetId(), ho_geometry);
1471 }
1472 MoFEMFunctionReturn(0);
1473}
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
a
constexpr double a
Definition: approx_sphere.cpp:46
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
PRESSURESET
@ PRESSURESET
Definition: definitions.h:165
FORCESET
@ FORCESET
Definition: definitions.h:164
NODESET
@ NODESET
Definition: definitions.h:159
SIDESET
@ SIDESET
Definition: definitions.h:160
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:49
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
lambda
constexpr double lambda
Definition: free_surface.cpp:101
get_f
auto get_f
Definition: free_surface.cpp:156
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
MoFEM::CoreInterface::add_ents_to_finite_element_by_dim
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:59
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:507
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:81
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:39
v
double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:52
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
ContactOps::normal
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:77
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:78
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
FTensor::dd
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
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MoFEM::Types::VectorDouble3
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:103
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEM::addHOOpsFace3D
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Definition: HODataOperators.hpp:517
MoFEM::invertTensor3by3
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1218
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1207
MoFEM::VectorDofs
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
Definition: EntitiesFieldData.hpp:35
ReactionDiffusionEquation::r
const double r
rate factor
Definition: reaction_diffusion.cpp:33
flux
double flux
impulse magnitude
Definition: photon_diffusion.cpp:64
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:31
_F
#define _F(n)
Definition: quad.c:25
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:95
MetaNeumannForces::setMassFluxOperators
static MoFEMErrorCode setMassFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Definition: SurfacePressure.cpp:1451
MetaNeumannForces::addNeumannFluxBCElements
static MoFEMErrorCode addNeumannFluxBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Definition: SurfacePressure.cpp:1426
MetaNeumannForces::addNeumannBCElements
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
Definition: SurfacePressure.cpp:1259
MetaNeumannForces::setMomentumFluxOperators
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
Definition: SurfacePressure.cpp:1354
MethodForForceScaling::applyScale
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Definition: MethodForForceScaling.hpp:33
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CubitMeshSets
this struct keeps basic methods for moab meshset about material and boundary conditions
Definition: BCMultiIndices.hpp:29
MoFEM::CubitMeshSets::getAttributes
MoFEMErrorCode getAttributes(std::vector< double > &attributes) const
get Cubit block attributes
Definition: BCMultiIndices.cpp:287
MoFEM::CubitMeshSets::getBcDataStructure
MoFEMErrorCode getBcDataStructure(CUBIT_BC_DATA_TYPE &data) const
Definition: BCMultiIndices.hpp:301
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:26
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:139
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:399
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1501
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1317
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1253
MoFEM::EntitiesFieldData::EntData::getFTensor1FieldData
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
Definition: EntitiesFieldData.hpp:1287
MoFEM::EntitiesFieldData::EntData::getFieldDofs
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
Definition: EntitiesFieldData.hpp:1268
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1213
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:95
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:114
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:103
NeumannForcesSurface::LinearVaringPresssure::getForce
MoFEMErrorCode getForce(const EntityHandle ent, const VectorDouble3 &coords, const VectorDouble3 &normal, VectorDouble3 &force)
Definition: SurfacePressure.cpp:26
NeumannForcesSurface::MyTriangleFE::MyTriangleFE
MyTriangleFE(MoFEM::Interface &m_field)
Definition: SurfacePressure.cpp:35
NeumannForcesSurface::OpCalculateDeformation
Operator for computing deformation gradients in side volumes.
Definition: SurfacePressure.hpp:320
NeumannForcesSurface::OpCalculateDeformation::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfacePressure.cpp:373
NeumannForcesSurface::OpGetTangent
Operator for computing tangent vectors.
Definition: SurfacePressure.hpp:236
NeumannForcesSurface::OpGetTangent::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfacePressure.cpp:236
NeumannForcesSurface::OpNeumannFlux::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfacePressure.cpp:938
NeumannForcesSurface::OpNeumannFlux::OpNeumannFlux
OpNeumannFlux(const std::string field_name, Vec _F, bCPressure &data, boost::ptr_vector< MethodForForceScaling > &methods_op, bool ho_geometry)
Definition: SurfacePressure.cpp:931
NeumannForcesSurface::OpNeumannForceAnalytical::OpNeumannForceAnalytical
OpNeumannForceAnalytical(const std::string field_name, Vec f, const Range tris, boost::ptr_vector< MethodForForceScaling > &methods_op, boost::shared_ptr< MethodForAnalyticalForce > &analytical_force_op, const bool ho_geometry=false)
Definition: SurfacePressure.cpp:108
NeumannForcesSurface::OpNeumannForceAnalytical::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfacePressure.cpp:118
NeumannForcesSurface::OpNeumannForce::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfacePressure.cpp:45
NeumannForcesSurface::OpNeumannForce::OpNeumannForce
OpNeumannForce(const std::string field_name, Vec _F, bCForce &data, boost::ptr_vector< MethodForForceScaling > &methods_op, bool ho_geometry=false)
Definition: SurfacePressure.cpp:38
NeumannForcesSurface::OpNeumannPressure
RHS-operator for pressure element (spatial configuration)
Definition: SurfacePressure.hpp:164
NeumannForcesSurface::OpNeumannPressure::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate pressure.
Definition: SurfacePressure.cpp:185
NeumannForcesSurface::OpNeumannPressure::OpNeumannPressure
OpNeumannPressure(const std::string field_name, Vec _F, bCPressure &data, boost::ptr_vector< MethodForForceScaling > &methods_op, bool ho_geometry=false)
Definition: SurfacePressure.cpp:178
NeumannForcesSurface::OpNeumannPressureLhs_dx_dX
LHS-operator for pressure element (spatial configuration)
Definition: SurfacePressure.hpp:256
NeumannForcesSurface::OpNeumannPressureLhs_dx_dX::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Compute left-hand side.
Definition: SurfacePressure.cpp:275
NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dX
LHS-operator for the pressure element (material configuration)
Definition: SurfacePressure.hpp:489
NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dX::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: SurfacePressure.cpp:528
NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Compute part of the left-hand side.
Definition: SurfacePressure.cpp:571
NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dx
LHS-operator for the pressure element (material configuration)
Definition: SurfacePressure.hpp:544
NeumannForcesSurface::OpNeumannPressureMaterialLhs_dX_dx::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: SurfacePressure.cpp:695
NeumannForcesSurface::OpNeumannPressureMaterialLhs::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Definition: SurfacePressure.cpp:647
NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX
RHS-operator for the pressure element (material configuration)
Definition: SurfacePressure.hpp:344
NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data)
Definition: SurfacePressure.cpp:440
NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data)
Definition: SurfacePressure.cpp:486
NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &row_data)
Integrate pressure in the material configuration.
Definition: SurfacePressure.cpp:409
NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs_dX_dX
LHS-operator (material configuration) on the side volume.
Definition: SurfacePressure.hpp:663
NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs_dX_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
Definition: SurfacePressure.cpp:864
NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs_dX_dx
LHS-operator (material configuration) on the side volume.
Definition: SurfacePressure.hpp:623
NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs_dX_dx::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
Definition: SurfacePressure.cpp:755
NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: SurfacePressure.cpp:717
NeumannForcesSurface::OpNeumannPressureMaterialVolOnSideLhs::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Definition: SurfacePressure.cpp:815
NeumannForcesSurface::mapPressure
std::map< int, bCPressure > mapPressure
Definition: SurfacePressure.hpp:107
NeumannForcesSurface::addPressureAle
MoFEMErrorCode addPressureAle(const std::string field_name_1, const std::string field_name_2, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, std::string side_fe_name, Vec F, Mat aij, int ms_id, bool ho_geometry=false, bool block_set=false)
Add operator to calculate pressure on element (in ALE)
Definition: SurfacePressure.cpp:1092
NeumannForcesSurface::mapForce
std::map< int, bCForce > mapForce
Definition: SurfacePressure.hpp:101
NeumannForcesSurface::addFlux
MoFEMErrorCode addFlux(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false)
Add flux element operator (integration on face)
Definition: SurfacePressure.cpp:1243
NeumannForcesSurface::methodsOp
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: SurfacePressure.hpp:109
NeumannForcesSurface::addPreassure
DEPRECATED MoFEMErrorCode addPreassure(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false, bool block_set=false)
Definition: SurfacePressure.cpp:1235
NeumannForcesSurface::addForce
MoFEMErrorCode addForce(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false, bool block_set=false)
Add operator to calculate forces on element.
Definition: SurfacePressure.cpp:987
NeumannForcesSurface::NeumannForcesSurface
NeumannForcesSurface(MoFEM::Interface &m_field)
Definition: SurfacePressure.hpp:93
NeumannForcesSurface::addLinearPressure
MoFEMErrorCode addLinearPressure(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false)
Add operator to calculate pressure on element.
Definition: SurfacePressure.cpp:1205
NeumannForcesSurface::addPressure
MoFEMErrorCode addPressure(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false, bool block_set=false)
Add operator to calculate pressure on element.
Definition: SurfacePressure.cpp:1049
NeumannForcesSurface::analyticalForceOp
boost::ptr_vector< MethodForAnalyticalForce > analyticalForceOp
Definition: SurfacePressure.hpp:110
NeumannForcesSurface::mField
MoFEM::Interface & mField
Definition: SurfacePressure.hpp:28
Generated on Tue Jul 5 2022 22:52:10 for MoFEM by Doxygen 1.9.4 and hosted at University of Glasgow