v0.14.0
FatPrismElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1 /** \file FatPrismElementForcesAndSourcesCore.cpp
2 
3 \brief Implementation of fat prism element
4 
5 */
6 
7 
8 
9 namespace MoFEM {
10 
11 FatPrismElementForcesAndSourcesCore::FatPrismElementForcesAndSourcesCore(
12  Interface &m_field)
13  : VolumeElementForcesAndSourcesCore(m_field, MBPRISM),
14  dataH1TrianglesOnly(MBPRISM), dataH1TroughThickness(MBPRISM),
15  opHOCoordsAndNormals(hoCoordsAtGaussPtsF3, nOrmals_at_GaussPtF3,
16  tAngent1_at_GaussPtF3, tAngent2_at_GaussPtF3,
17  hoCoordsAtGaussPtsF4, nOrmals_at_GaussPtF4,
18  tAngent1_at_GaussPtF4, tAngent2_at_GaussPtF4) {
19  CHK_THROW_MESSAGE(createDataOnElement(MBPRISM),
20  "Problem with creation data on element");
21 }
22 
23 MoFEMErrorCode FatPrismElementForcesAndSourcesCore::operator()() {
24  MoFEMFunctionBegin;
25 
26  if (numeredEntFiniteElementPtr->getEntType() != MBPRISM)
27  MoFEMFunctionReturnHot(0);
28 
29  auto get_fe_coordinates = [&]() {
30  MoFEMFunctionBegin;
31  EntityHandle ent = numeredEntFiniteElementPtr->getEnt();
32  int num_nodes;
33  const EntityHandle *conn;
34  CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
35  coords.resize(num_nodes * 3, false);
36  CHKERR mField.get_moab().get_coords(conn, num_nodes,
37  &*coords.data().begin());
38  MoFEMFunctionReturn(0);
39  };
40 
41  auto calculate_area_of_triangles = [&] {
42  MoFEMFunctionBegin;
43  normal.resize(6, false);
44  CHKERR Tools::getTriNormal(&coords[0], &normal[0]);
45  CHKERR Tools::getTriNormal(&coords[9], &normal[3]);
46  aRea[0] = cblas_dnrm2(3, &normal[0], 1) * 0.5;
47  aRea[1] = cblas_dnrm2(3, &normal[3], 1) * 0.5;
48  MoFEMFunctionReturn(0);
49  };
50 
51  CHKERR get_fe_coordinates();
52  CHKERR calculate_area_of_triangles();
53 
54  CHKERR getSpacesAndBaseOnEntities(dataH1);
55  CHKERR getSpacesAndBaseOnEntities(dataH1TrianglesOnly);
56  CHKERR getSpacesAndBaseOnEntities(dataH1TroughThickness);
57 
58  auto get_h1_base_data = [&](auto &dataH1) {
59  MoFEMFunctionBegin;
60  CHKERR getEntitySense<MBEDGE>(dataH1);
61  CHKERR getEntitySense<MBTRI>(dataH1);
62  CHKERR getEntitySense<MBQUAD>(dataH1);
63  CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
64  CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
65  CHKERR getEntityDataOrder<MBQUAD>(dataH1, H1);
66  CHKERR getEntityDataOrder<MBPRISM>(dataH1, H1);
67  // Triangles only
68  CHKERR getEntitySense<MBEDGE>(dataH1TrianglesOnly);
69  CHKERR getEntitySense<MBTRI>(dataH1TrianglesOnly);
70  CHKERR getEntityDataOrder<MBEDGE>(dataH1TrianglesOnly, H1);
71  CHKERR getEntityDataOrder<MBTRI>(dataH1TrianglesOnly, H1);
72  // Through thickness
73  CHKERR getEntitySense<MBEDGE>(dataH1TroughThickness);
74  CHKERR getEntityDataOrder<MBEDGE>(dataH1TroughThickness, H1);
75  MoFEMFunctionReturn(0);
76  };
77 
78  // H1
79  if ((dataH1.spacesOnEntities[MBEDGE]).test(H1))
80  CHKERR get_h1_base_data(dataH1);
81 
82  // Hdiv
83  if ((dataH1.spacesOnEntities[MBTRI]).test(HDIV))
84  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented yet");
85 
86  // Hcurl
87  if ((dataH1.spacesOnEntities[MBEDGE]).test(HCURL))
88  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented yet");
89 
90  // L2
91  if ((dataH1.spacesOnEntities[MBPRISM]).test(L2))
92  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented yet");
93 
94  // get approx. on triangles, i.e. faces 3 and 4
95  auto set_gauss_points_on_triangle = [&](int &nb_gauss_pts_on_faces) {
96  MoFEMFunctionBegin;
97  int order_triangles_only = 1;
98  int valid_edges[] = {1, 1, 1, 0, 0, 0, 1, 1, 1};
99  for (unsigned int ee = 0; ee < 9; ee++) {
100  if (!valid_edges[ee])
101  continue;
102  order_triangles_only = std::max(
103  order_triangles_only,
104  dataH1TrianglesOnly.dataOnEntities[MBEDGE][ee].getOrder());
105  }
106  for (unsigned int ff = 3; ff <= 4; ff++) {
107  order_triangles_only = std::max(
108  order_triangles_only,
109  dataH1TrianglesOnly.dataOnEntities[MBTRI][ff].getOrder());
110  }
111  for (unsigned int qq = 0; qq < 3; qq++) {
112  order_triangles_only = std::max(
113  order_triangles_only,
114  dataH1TroughThickness.dataOnEntities[MBQUAD][qq].getOrder());
115  }
116  order_triangles_only = std::max(
117  order_triangles_only,
118  dataH1TroughThickness.dataOnEntities[MBPRISM][0].getOrder());
119 
120  // integration pts on the triangles surfaces
121  nb_gauss_pts_on_faces = 0;
122  int rule = getRuleTrianglesOnly(order_triangles_only);
123  if (rule >= 0) {
124  if (rule < QUAD_2D_TABLE_SIZE) {
125  if (QUAD_2D_TABLE[rule]->dim != 2) {
126  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
127  }
128  if (QUAD_2D_TABLE[rule]->order < rule) {
129  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
130  "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
131  }
132  nb_gauss_pts_on_faces = QUAD_2D_TABLE[rule]->npoints;
133  gaussPtsTrianglesOnly.resize(3, nb_gauss_pts_on_faces, false);
134  cblas_dcopy(nb_gauss_pts_on_faces, &QUAD_2D_TABLE[rule]->points[1], 3,
135  &gaussPtsTrianglesOnly(0, 0), 1);
136  cblas_dcopy(nb_gauss_pts_on_faces, &QUAD_2D_TABLE[rule]->points[2], 3,
137  &gaussPtsTrianglesOnly(1, 0), 1);
138  cblas_dcopy(nb_gauss_pts_on_faces, QUAD_2D_TABLE[rule]->weights, 1,
139  &gaussPtsTrianglesOnly(2, 0), 1);
140  dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(
141  nb_gauss_pts_on_faces, 3, false);
142  double *shape_ptr = &*dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0]
143  .getN(NOBASE)
144  .data()
145  .begin();
146  cblas_dcopy(3 * nb_gauss_pts_on_faces, QUAD_2D_TABLE[rule]->points, 1,
147  shape_ptr, 1);
148  dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(
149  1, 6, false);
150  std::copy(Tools::diffShapeFunMBTRI.begin(),
151  Tools::diffShapeFunMBTRI.end(),
152  &*dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0]
153  .getDiffN(NOBASE)
154  .data()
155  .begin());
156  } else
157  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
158  "rule > quadrature order %d < %d", rule, QUAD_2D_TABLE_SIZE);
159 
160  } else {
161  CHKERR setGaussPtsTrianglesOnly(order_triangles_only);
162  nb_gauss_pts_on_faces = gaussPtsTrianglesOnly.size2();
163  if (nb_gauss_pts_on_faces == 0)
164  MoFEMFunctionReturnHot(0);
165  dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(
166  nb_gauss_pts_on_faces, 3, false);
167  if (nb_gauss_pts_on_faces) {
168  CHKERR ShapeMBTRI(&*dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0]
169  .getN(NOBASE)
170  .data()
171  .begin(),
172  &gaussPtsTrianglesOnly(0, 0),
173  &gaussPtsTrianglesOnly(1, 0), nb_gauss_pts_on_faces);
174  dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(
175  1, 6, false);
176  std::copy(Tools::diffShapeFunMBTRI.begin(),
177  Tools::diffShapeFunMBTRI.end(),
178  &*dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0]
179  .getDiffN(NOBASE)
180  .data()
181  .begin());
182  }
183  }
184  MoFEMFunctionReturn(0);
185  };
186 
187  // approx. trough prism thickness
188  auto set_gauss_points_through_thickness =
189  [&](int &nb_gauss_pts_through_thickness) {
190  MoFEMFunctionBegin;
191  nb_gauss_pts_through_thickness = 0;
192  int order_thickness = 1;
193  for (unsigned int ee = 3; ee <= 5; ee++) {
194  order_thickness = std::max(
195  order_thickness,
196  dataH1TroughThickness.dataOnEntities[MBEDGE][ee].getOrder());
197  }
198  for (unsigned int qq = 0; qq < 3; qq++) {
199  order_thickness = std::max(
200  order_thickness,
201  dataH1TroughThickness.dataOnEntities[MBQUAD][qq].getOrder());
202  }
203  order_thickness = std::max(
204  order_thickness,
205  dataH1TroughThickness.dataOnEntities[MBPRISM][0].getOrder());
206  // integration points
207  int rule = getRuleThroughThickness(order_thickness);
208  if (rule >= 0) {
209  if (rule < QUAD_1D_TABLE_SIZE) {
210  if (QUAD_1D_TABLE[rule]->dim != 1) {
211  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
212  "wrong dimension");
213  }
214  if (QUAD_1D_TABLE[rule]->order < rule) {
215  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
216  "wrong order %d != %d", QUAD_1D_TABLE[rule]->order,
217  rule);
218  }
219  nb_gauss_pts_through_thickness = QUAD_1D_TABLE[rule]->npoints;
220  gaussPtsThroughThickness.resize(2, nb_gauss_pts_through_thickness,
221  false);
222  cblas_dcopy(nb_gauss_pts_through_thickness,
223  &QUAD_1D_TABLE[rule]->points[1], 2,
224  &gaussPtsThroughThickness(0, 0), 1);
225  cblas_dcopy(nb_gauss_pts_through_thickness,
226  QUAD_1D_TABLE[rule]->weights, 1,
227  &gaussPtsThroughThickness(1, 0), 1);
228  } else {
229  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
230  "rule > quadrature order %d < %d", rule,
231  QUAD_1D_TABLE_SIZE);
232  nb_gauss_pts_through_thickness = 0;
233  }
234  } else {
235  CHKERR setGaussPtsThroughThickness(order_thickness);
236  nb_gauss_pts_through_thickness = gaussPtsThroughThickness.size2();
237  }
238  MoFEMFunctionReturn(0);
239  };
240 
241  // Generate integration pts.
242  auto set_gauss_points_in_volume = [&](int nb_gauss_pts_on_faces,
243  int nb_gauss_pts_through_thickness,
244  int &nb_gauss_pts) {
245  MoFEMFunctionBegin;
246  nb_gauss_pts = nb_gauss_pts_on_faces * nb_gauss_pts_through_thickness;
247  gaussPts.resize(4, nb_gauss_pts, false);
248  int gg = 0;
249  for (int ggf = 0; ggf < nb_gauss_pts_on_faces; ggf++) {
250  for (int ggt = 0; ggt < nb_gauss_pts_through_thickness; ggt++, gg++) {
251  gaussPts(0, gg) = gaussPtsTrianglesOnly(0, ggf);
252  gaussPts(1, gg) = gaussPtsTrianglesOnly(1, ggf);
253  gaussPts(2, gg) = gaussPtsThroughThickness(0, ggt);
254  gaussPts(3, gg) =
255  gaussPtsTrianglesOnly(2, ggf) * gaussPtsThroughThickness(1, ggt);
256  }
257  }
258  MoFEMFunctionReturn(0);
259  };
260 
261  int nb_gauss_pts, nb_gauss_pts_through_thickness, nb_gauss_pts_on_faces;
262  CHKERR set_gauss_points_on_triangle(nb_gauss_pts_on_faces);
263  if (!nb_gauss_pts_on_faces)
264  MoFEMFunctionReturnHot(0);
265  CHKERR set_gauss_points_through_thickness(nb_gauss_pts_through_thickness);
266  if (!nb_gauss_pts_through_thickness)
267  MoFEMFunctionReturnHot(0);
268  CHKERR set_gauss_points_in_volume(
269  nb_gauss_pts_on_faces, nb_gauss_pts_through_thickness, nb_gauss_pts);
270 
271  auto calc_coordinates_at_triangles = [&]() {
272  MoFEMFunctionBegin;
273  coordsAtGaussPtsTrianglesOnly.resize(nb_gauss_pts_on_faces, 6, false);
274  for (int gg = 0; gg < nb_gauss_pts_on_faces; gg++) {
275  for (int dd = 0; dd < 3; dd++) {
276  coordsAtGaussPtsTrianglesOnly(gg, dd) =
277  cblas_ddot(3,
278  &dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getN(
279  NOBASE)(gg, 0),
280  1, &coords[dd], 3);
281  coordsAtGaussPtsTrianglesOnly(gg, 3 + dd) =
282  cblas_ddot(3,
283  &dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getN(
284  NOBASE)(gg, 0),
285  1, &coords[9 + dd], 3);
286  }
287  }
288  MoFEMFunctionReturn(0);
289  };
290 
291  auto calc_vertex_base_on_prism = [&]() {
292  MoFEMFunctionBegin;
293  // Calculate "nobase" base functions on prism, this is cartesian product
294  // of base functions on triangles with base functions through thickness
295  // FIXME: This could be effectively implemented with tensors
296  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 6,
297  false);
298  dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(nb_gauss_pts, 18,
299  false);
300  for (int dd = 0; dd != 6; dd++) {
301  int gg = 0;
302  for (int ggf = 0; ggf < nb_gauss_pts_on_faces; ggf++) {
303  int ddd = dd > 2 ? dd - 3 : dd;
304  double tri_n = dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getN(
305  NOBASE)(ggf, ddd);
306  double dksi_tri_n =
307  dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(
308  0, 2 * ddd + 0);
309  double deta_tri_n =
310  dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(
311  0, 2 * ddd + 1);
312  for (int ggt = 0; ggt < nb_gauss_pts_through_thickness; ggt++, gg++) {
313  double zeta = gaussPtsThroughThickness(0, ggt);
314  double dzeta, edge_shape;
315  if (dd < 3) {
316  dzeta = diffN_MBEDGE0;
317  edge_shape = N_MBEDGE0(zeta);
318  } else {
319  dzeta = diffN_MBEDGE1;
320  edge_shape = N_MBEDGE1(zeta);
321  }
322  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, dd) =
323  tri_n * edge_shape;
324  dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(gg, 3 * dd + 0) =
325  dksi_tri_n * edge_shape;
326  dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(gg, 3 * dd + 1) =
327  deta_tri_n * edge_shape;
328  dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(gg, 3 * dd + 2) =
329  tri_n * dzeta;
330  }
331  }
332  }
333  MoFEMFunctionReturn(0);
334  };
335 
336  auto calc_base_on_prism = [&]() {
337  MoFEMFunctionBegin;
338  // Calculate base functions on prism
339  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
340  if (dataH1.bAse.test(b)) {
341  switch (static_cast<FieldApproximationBase>(b)) {
342  case AINSWORTH_LEGENDRE_BASE:
343  case AINSWORTH_LOBATTO_BASE:
344  if (dataH1.spacesOnEntities[MBVERTEX].test(H1)) {
345  CHKERR FatPrismPolynomialBase().getValue(
346  gaussPts,
347  boost::shared_ptr<BaseFunctionCtx>(
348  new FatPrismPolynomialBaseCtx(
349  dataH1, dataH1TrianglesOnly, dataH1TroughThickness,
350  gaussPtsTrianglesOnly, gaussPtsThroughThickness,
351  mField.get_moab(), numeredEntFiniteElementPtr.get(), H1,
352  static_cast<FieldApproximationBase>(b), NOBASE)));
353  }
354  if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
355  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
356  "Not yet implemented");
357  }
358  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
359  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
360  "Not yet implemented");
361  }
362  if (dataH1.spacesOnEntities[MBTET].test(L2)) {
363  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
364  "Not yet implemented");
365  }
366  break;
367  default:
368  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
369  "Not yet implemented");
370  }
371  }
372  }
373  MoFEMFunctionReturn(0);
374  };
375 
376  auto calc_coordinate_on_prism = [&]() {
377  MoFEMFunctionBegin;
378  /// Calculate coordinates at integration points
379  coordsAtGaussPts.resize(nb_gauss_pts, 3, false);
380  for (int gg = 0; gg < nb_gauss_pts; gg++) {
381  for (int dd = 0; dd < 3; dd++) {
382  coordsAtGaussPts(gg, dd) = cblas_ddot(
383  6, &dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, 0), 1,
384  &coords[dd], 3);
385  }
386  }
387  MoFEMFunctionReturn(0);
388  };
389 
390  auto calculate_volume = [&]() {
391  auto get_t_w = [&] {
392  return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(&gaussPts(2, 0));
393  };
394 
395  auto get_t_coords = [&]() {
396  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
397  &coords[0], &coords[1], &coords[2]);
398  };
399 
400  FTensor::Index<'i', 3> i;
401  FTensor::Index<'j', 3> j;
402  FTensor::Tensor2<double, 3, 3> t_jac;
403 
404  const size_t nb_gauss_pts = gaussPts.size2();
405  auto t_diff_n =
406  dataH1.dataOnEntities[MBVERTEX][0].getFTensor1DiffN<3>(NOBASE);
407 
408  double vol = 0;
409  auto t_w = get_t_w();
410  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
411 
412  auto t_coords = get_t_coords();
413  t_jac(i, j) = 0;
414  for (size_t n = 0; n != 6; ++n) {
415  t_jac(i, j) += t_coords(i) * t_diff_n(j);
416  ++t_diff_n;
417  ++t_coords;
418  }
419 
420  double det;
421  CHKERR determinantTensor3by3(t_jac, det);
422  vol += det * t_w / 2;
423 
424  ++t_w;
425  }
426 
427  return vol;
428  };
429 
430  auto calc_ho_triangle_face_normals = [&]() {
431  MoFEMFunctionBegin;
432 
433  auto check_field = [&]() {
434  auto field_it =
435  fieldsPtr->get<FieldName_mi_tag>().find(meshPositionsFieldName);
436  if (field_it != fieldsPtr->get<FieldName_mi_tag>().end())
437  if ((numeredEntFiniteElementPtr->getBitFieldIdData() &
438  (*field_it)->getId())
439  .any())
440  return true;
441  return false;
442  };
443 
444  // Check if field meshPositionsFieldName exist
445  if (check_field()) {
446  hoCoordsAtGaussPtsF3.resize(nb_gauss_pts_on_faces, 3, false);
447  nOrmals_at_GaussPtF3.resize(nb_gauss_pts_on_faces, 3, false);
448  tAngent1_at_GaussPtF3.resize(nb_gauss_pts_on_faces, 3, false);
449  tAngent2_at_GaussPtF3.resize(nb_gauss_pts_on_faces, 3, false);
450  hoCoordsAtGaussPtsF4.resize(nb_gauss_pts_on_faces, 3, false);
451  nOrmals_at_GaussPtF4.resize(nb_gauss_pts_on_faces, 3, false);
452  tAngent1_at_GaussPtF4.resize(nb_gauss_pts_on_faces, 3, false);
453  tAngent2_at_GaussPtF4.resize(nb_gauss_pts_on_faces, 3, false);
454  const auto bit_number =
455  mField.get_field_bit_number(meshPositionsFieldName);
456  CHKERR getNodesFieldData(dataH1TrianglesOnly, bit_number);
457  CHKERR getEntityFieldData(dataH1TrianglesOnly, bit_number, MBEDGE);
458  CHKERR getEntityFieldData(dataH1TrianglesOnly, bit_number, MBEDGE);
459  CHKERR opHOCoordsAndNormals.opRhs(dataH1TrianglesOnly);
460  CHKERR opHOCoordsAndNormals.calculateNormals();
461  } else {
462  hoCoordsAtGaussPtsF3.resize(0, 0, false);
463  nOrmals_at_GaussPtF3.resize(0, 0, false);
464  tAngent1_at_GaussPtF3.resize(0, 0, false);
465  tAngent2_at_GaussPtF3.resize(0, 0, false);
466  hoCoordsAtGaussPtsF4.resize(0, 0, false);
467  nOrmals_at_GaussPtF4.resize(0, 0, false);
468  tAngent1_at_GaussPtF4.resize(0, 0, false);
469  tAngent2_at_GaussPtF4.resize(0, 0, false);
470  }
471  MoFEMFunctionReturn(0);
472  };
473 
474  CHKERR calc_coordinates_at_triangles();
475  CHKERR calc_vertex_base_on_prism();
476  CHKERR calc_base_on_prism();
477  CHKERR calc_coordinate_on_prism();
478  vOlume = calculate_volume();
479  CHKERR calc_ho_triangle_face_normals();
480 
481  // Iterate over operators
482  CHKERR loopOverOperators();
483 
484  MoFEMFunctionReturn(0);
485 }
486 
487 MoFEMErrorCode FatPrismElementForcesAndSourcesCore::UserDataOperator::setPtrFE(
488  ForcesAndSourcesCore *ptr) {
489  MoFEMFunctionBeginHot;
490  if (!(ptrFE = dynamic_cast<FatPrismElementForcesAndSourcesCore *>(ptr)))
491  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
492  "User operator and finite element do not work together");
493  MoFEMFunctionReturnHot(0);
494 }
495 
496 } // namespace MoFEM
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
diffN_MBEDGE0
#define diffN_MBEDGE0
derivative of edge shape function
Definition: fem_tools.h:107
MoFEM::FatPrismElementForcesAndSourcesCore::dataH1TrianglesOnly
EntitiesFieldData dataH1TrianglesOnly
Definition: FatPrismElementForcesAndSourcesCore.hpp:193
MoFEM::FatPrismElementForcesAndSourcesCore::hoCoordsAtGaussPtsF3
MatrixDouble hoCoordsAtGaussPtsF3
Definition: FatPrismElementForcesAndSourcesCore.hpp:196
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::VolumeElementForcesAndSourcesCore::meshPositionsFieldName
std::string meshPositionsFieldName
Definition: VolumeElementForcesAndSourcesCore.hpp:32
MoFEM::DataOperator::opRhs
virtual MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
Definition: DataOperators.cpp:140
MoFEM::FatPrismPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: FatPrismPolynomialBase.cpp:50
LASTBASE
@ LASTBASE
Definition: definitions.h:69
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
MoFEM::FatPrismElementForcesAndSourcesCore::hoCoordsAtGaussPtsF4
MatrixDouble hoCoordsAtGaussPtsF4
Definition: FatPrismElementForcesAndSourcesCore.hpp:200
NOBASE
@ NOBASE
Definition: definitions.h:59
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::FatPrismElementForcesAndSourcesCore::tAngent1_at_GaussPtF4
MatrixDouble tAngent1_at_GaussPtF4
Definition: FatPrismElementForcesAndSourcesCore.hpp:202
MoFEM::FatPrismElementForcesAndSourcesCore::tAngent1_at_GaussPtF3
MatrixDouble tAngent1_at_GaussPtF3
Definition: FatPrismElementForcesAndSourcesCore.hpp:198
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
QUAD_2D_TABLE
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
MoFEM::FatPrismElementForcesAndSourcesCore::aRea
double aRea[2]
Definition: FatPrismElementForcesAndSourcesCore.hpp:186
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::FatPrismElementForcesAndSourcesCore::tAngent2_at_GaussPtF4
MatrixDouble tAngent2_at_GaussPtF4
Definition: FatPrismElementForcesAndSourcesCore.hpp:203
MoFEM::FatPrismElementForcesAndSourcesCore::gaussPtsTrianglesOnly
MatrixDouble gaussPtsTrianglesOnly
Definition: FatPrismElementForcesAndSourcesCore.hpp:189
MoFEM::FieldName_mi_tag
MultiIndex Tag for field name.
Definition: TagMultiIndices.hpp:67
N_MBEDGE0
#define N_MBEDGE0(x)
edge shape function
Definition: fem_tools.h:105
QUAD_2D_TABLE_SIZE
#define QUAD_2D_TABLE_SIZE
Definition: quad.h:174
MoFEM::FatPrismPolynomialBase
Calculate base functions on tetrahedral.
Definition: FatPrismPolynomialBase.hpp:56
MoFEM::EntitiesFieldData::spacesOnEntities
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Definition: EntitiesFieldData.hpp:49
MoFEM::FatPrismElementForcesAndSourcesCore::coordsAtGaussPtsTrianglesOnly
MatrixDouble coordsAtGaussPtsTrianglesOnly
Definition: FatPrismElementForcesAndSourcesCore.hpp:190
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:126
FTensor::Tensor2< double, 3, 3 >
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::ForcesAndSourcesCore::getEntityFieldData
MoFEMErrorCode getEntityFieldData(EntitiesFieldData &data, const int bit_number, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
Definition: ForcesAndSourcesCore.cpp:770
MoFEM::VolumeElementForcesAndSourcesCore::coords
VectorDouble coords
Definition: VolumeElementForcesAndSourcesCore.hpp:89
MoFEM::BasicMethod::fieldsPtr
const Field_multiIndex * fieldsPtr
raw pointer to fields container
Definition: LoopMethods.hpp:252
QUAD_1D_TABLE
static QUAD *const QUAD_1D_TABLE[]
Definition: quad.h:164
MoFEM::FatPrismElementForcesAndSourcesCore::setGaussPtsThroughThickness
virtual MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
Definition: FatPrismElementForcesAndSourcesCore.hpp:45
MoFEM::FatPrismElementForcesAndSourcesCore::FatPrismElementForcesAndSourcesCore
FatPrismElementForcesAndSourcesCore(Interface &m_field)
Definition: FatPrismElementForcesAndSourcesCore.cpp:11
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::FatPrismElementForcesAndSourcesCore::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: FatPrismElementForcesAndSourcesCore.cpp:23
MoFEM::VolumeElementForcesAndSourcesCore::num_nodes
int num_nodes
Definition: VolumeElementForcesAndSourcesCore.hpp:100
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
N_MBEDGE1
#define N_MBEDGE1(x)
edge shape function
Definition: fem_tools.h:106
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::ForcesAndSourcesCore::coordsAtGaussPts
MatrixDouble coordsAtGaussPts
coordinated at gauss points
Definition: ForcesAndSourcesCore.hpp:544
MoFEM::FatPrismElementForcesAndSourcesCore::dataH1TroughThickness
EntitiesFieldData dataH1TroughThickness
Definition: FatPrismElementForcesAndSourcesCore.hpp:194
MoFEM::VolumeElementForcesAndSourcesCore::conn
const EntityHandle * conn
Definition: VolumeElementForcesAndSourcesCore.hpp:101
MoFEM::ForcesAndSourcesCore::createDataOnElement
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
Definition: ForcesAndSourcesCore.cpp:1337
MoFEM::FatPrismPolynomialBaseCtx
Class used to pass element data to calculate base functions on fat prism.
Definition: FatPrismPolynomialBase.hpp:23
QUAD_::npoints
int npoints
Definition: quad.h:29
MoFEM::FatPrismElementForcesAndSourcesCore::getRuleTrianglesOnly
virtual int getRuleTrianglesOnly(int order)
Definition: FatPrismElementForcesAndSourcesCore.hpp:36
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
MoFEM::Tools::diffShapeFunMBTRI
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:104
MoFEM::FatPrismElementForcesAndSourcesCore::setGaussPtsTrianglesOnly
virtual MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
Definition: FatPrismElementForcesAndSourcesCore.hpp:39
MoFEM::FatPrismElementForcesAndSourcesCore
FatPrism finite element.
Definition: FatPrismElementForcesAndSourcesCore.hpp:31
MoFEM::EntitiesFieldData::bAse
std::bitset< LASTBASE > bAse
bases on element
Definition: EntitiesFieldData.hpp:42
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::ForcesAndSourcesCore::getSpacesAndBaseOnEntities
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
Definition: ForcesAndSourcesCore.cpp:990
FTensor::Index< 'i', 3 >
diffN_MBEDGE1
#define diffN_MBEDGE1
derivative of edge shape function
Definition: fem_tools.h:108
convert.n
n
Definition: convert.py:82
MoFEM::FatPrismElementForcesAndSourcesCore::gaussPtsThroughThickness
MatrixDouble gaussPtsThroughThickness
Definition: FatPrismElementForcesAndSourcesCore.hpp:191
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
MoFEM::VolumeElementForcesAndSourcesCore::vOlume
double & vOlume
Definition: VolumeElementForcesAndSourcesCore.hpp:98
MoFEM::ForcesAndSourcesCore::dataH1
EntitiesFieldData & dataH1
Definition: ForcesAndSourcesCore.hpp:471
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
QUAD_1D_TABLE_SIZE
#define QUAD_1D_TABLE_SIZE
Definition: quad.h:163
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
FTensor::Tensor0
Definition: Tensor0.hpp:16
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::FatPrismElementForcesAndSourcesCore::getRuleThroughThickness
virtual int getRuleThroughThickness(int order)
Definition: FatPrismElementForcesAndSourcesCore.hpp:37
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::ForcesAndSourcesCore::gaussPts
MatrixDouble gaussPts
Matrix of integration points.
Definition: ForcesAndSourcesCore.hpp:109
MoFEM::FatPrismElementForcesAndSourcesCore::nOrmals_at_GaussPtF4
MatrixDouble nOrmals_at_GaussPtF4
Definition: FatPrismElementForcesAndSourcesCore.hpp:201
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: EntitiesFieldData.hpp:57
MoFEM::ForcesAndSourcesCore::loopOverOperators
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
Definition: ForcesAndSourcesCore.cpp:1379
MoFEM::FEMethod::numeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Definition: LoopMethods.hpp:383
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEM::field_it
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it
Definition: Projection10NodeCoordsOnField.cpp:123
MoFEM::ForcesAndSourcesCore::getNodesFieldData
MoFEMErrorCode getNodesFieldData(EntitiesFieldData &data, const int bit_number) const
Get data on nodes.
Definition: ForcesAndSourcesCore.cpp:642
MoFEM::FatPrismElementForcesAndSourcesCore::tAngent2_at_GaussPtF3
MatrixDouble tAngent2_at_GaussPtF3
Definition: FatPrismElementForcesAndSourcesCore.hpp:199
MoFEM::FatPrismElementForcesAndSourcesCore::nOrmals_at_GaussPtF3
MatrixDouble nOrmals_at_GaussPtF3
Definition: FatPrismElementForcesAndSourcesCore.hpp:197
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MoFEM::Tools::getTriNormal
static MoFEMErrorCode getTriNormal(const double *coords, double *normal, double *d_normal=nullptr)
Get the Tri Normal objectGet triangle normal.
Definition: Tools.cpp:353
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::FatPrismElementForcesAndSourcesCore::normal
VectorDouble normal
Definition: FatPrismElementForcesAndSourcesCore.hpp:187
MoFEM::FatPrismElementForcesAndSourcesCore::opHOCoordsAndNormals
OpGetCoordsAndNormalsOnPrism opHOCoordsAndNormals
Definition: FatPrismElementForcesAndSourcesCore.hpp:204
MoFEM::ForcesAndSourcesCore::UserDataOperator::ptrFE
ForcesAndSourcesCore * ptrFE
Definition: ForcesAndSourcesCore.hpp:985
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::ForcesAndSourcesCore::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:24
MoFEM::OpGetCoordsAndNormalsOnPrism::calculateNormals
MoFEMErrorCode calculateNormals()
Definition: DataOperators.cpp:472
ShapeMBTRI
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:182
MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator::setPtrFE
MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr)
Definition: FatPrismElementForcesAndSourcesCore.cpp:487
Generated by Doxygen 1.8.17 and hosted at University of Glasgow