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