v0.9.0
FaceElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1 /** \file FaceElementForcesAndSourcesCore.cpp
2 
3 \brief Implementation of face element
4 
5 */
6 
7 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 namespace MoFEM {
22 
24  Interface &m_field)
25  : ForcesAndSourcesCore(m_field),
26  meshPositionsFieldName("MESH_NODE_POSITIONS"),
27  opHOCoordsAndNormals(hoCoordsAtGaussPts, normalsAtGaussPts,
28  tangentOneAtGaussPts, tangentTwoAtGaussPts),
29  opContravariantTransform(nOrmal, normalsAtGaussPts),
30  opCovariantTransform(nOrmal, normalsAtGaussPts, tangentOne,
31  tangentOneAtGaussPts, tangentTwo,
32  tangentTwoAtGaussPts) {
33 }
34 
38 
39  auto type = numeredEntFiniteElementPtr->getEntType();
40 
41  if (type == MBQUAD) {
42 
44  CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
45  coords.resize(num_nodes * 3, false);
46  CHKERR mField.get_moab().get_coords(conn, num_nodes,
47  &*coords.data().begin());
48 
49  const size_t nb_gauss_pts = gaussPts.size2();
50  normalsAtGaussPts.resize(nb_gauss_pts, 3);
51  tangentOneAtGaussPts.resize(nb_gauss_pts, 3);
52  tangentTwoAtGaussPts.resize(nb_gauss_pts, 3);
53  normalsAtGaussPts.clear();
54  tangentOneAtGaussPts.clear();
55  tangentTwoAtGaussPts.clear();
56 
57  auto get_ftensor_from_mat_3d = [](MatrixDouble &m) {
59  &m(0, 0), &m(0, 1), &m(0, 2));
60  };
61  auto get_ftensor_from_mat_2d = [](MatrixDouble &m) {
63  &m(0, 1));
64  };
65  auto get_ftensor_from_vec_3d = [](VectorDouble &v) {
66  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
67  &v[2]);
68  };
69 
70  auto t_t1 = get_ftensor_from_mat_3d(tangentOneAtGaussPts);
71  auto t_t2 = get_ftensor_from_mat_3d(tangentTwoAtGaussPts);
72  auto t_normal = get_ftensor_from_mat_3d(normalsAtGaussPts);
73 
74  auto t_diff = get_ftensor_from_mat_2d(
75  dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE));
76 
80 
83 
84  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
85  auto t_coords = get_ftensor_from_vec_3d(coords);
86 
87  for (int nn = 0; nn != num_nodes; ++nn) {
88  t_t1(i) += t_coords(i) * t_diff(N0);
89  t_t2(i) += t_coords(i) * t_diff(N1);
90  ++t_diff;
91  ++t_coords;
92  }
93  t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
94 
95  ++t_t1;
96  ++t_t2;
97  ++t_normal;
98  }
99  }
100 
102 }
103 
106 
108  CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
109  coords.resize(num_nodes * 3, false);
110  CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
111  nOrmal.resize(3, false);
112  tangentOne.resize(3, false);
113  tangentTwo.resize(3, false);
114 
115  auto calc_normal = [&](const double *diff_ptr) {
118  &coords[0], &coords[1], &coords[2]);
120  &nOrmal[0], &nOrmal[1], &nOrmal[2]);
122  &tangentOne[0], &tangentOne[1], &tangentOne[2]);
124  &tangentTwo[0], &tangentTwo[1], &tangentTwo[2]);
126  diff_ptr, &diff_ptr[1]);
127 
131 
134  t_t1(i) = 0;
135  t_t2(i) = 0;
136 
137  for (int nn = 0; nn != num_nodes; ++nn) {
138  t_t1(i) += t_coords(i) * t_diff(N0);
139  t_t2(i) += t_coords(i) * t_diff(N1);
140  ++t_coords;
141  ++t_diff;
142  }
143  t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
144  aRea = sqrt(t_normal(i) * t_normal(i));
146  };
147 
148  const double *diff_ptr;
149  switch (numeredEntFiniteElementPtr->getEntType()) {
150  case MBTRI:
151  diff_ptr = Tools::diffShapeFunMBTRI.data();
152  CHKERR calc_normal(diff_ptr);
153  //FIXME: Normal should be divided not the area for triangle!!
154  aRea /= 2;
155  break;
156  case MBQUAD:
157  diff_ptr = Tools::diffShapeFunMBQUADAtCenter.data();
158  CHKERR calc_normal(diff_ptr);
159  break;
160  default:
161  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Element type not implemented");
162  }
163 
165 }
166 
169  // Set integration points
170  int order_data = getMaxDataOrder();
171  int order_row = getMaxRowOrder();
172  int order_col = getMaxColOrder();
173  int rule = getRule(order_row, order_col, order_data);
174 
175  auto set_integration_pts_for_tri = [&]() {
177  if (rule < QUAD_2D_TABLE_SIZE) {
178  if (QUAD_2D_TABLE[rule]->dim != 2) {
179  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
180  }
181  if (QUAD_2D_TABLE[rule]->order < rule) {
182  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
183  "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
184  }
185  const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
186  gaussPts.resize(3, nb_gauss_pts, false);
187  cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
188  &gaussPts(0, 0), 1);
189  cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
190  &gaussPts(1, 0), 1);
191  cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
192  &gaussPts(2, 0), 1);
193  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
194  false);
195  double *shape_ptr =
196  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
197  cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
198  1);
199  dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(3, 2, false);
200  std::copy(
202  dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
203  }
204  else {
205  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
206  "rule > quadrature order %d < %d", rule, QUAD_2D_TABLE_SIZE);
207  }
209  };
210 
211  auto calc_base_for_quad = [&]() {
213  const size_t nb_gauss_pts = gaussPts.size2();
214  auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
215  auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
216  base.resize(nb_gauss_pts, 4, false);
217  diff_base.resize(nb_gauss_pts, 8, false);
218  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
219  const double ksi = gaussPts(0, gg);
220  const double zeta = gaussPts(1, gg);
221  base(gg, 0) = N_MBQUAD0(ksi, zeta);
222  base(gg, 1) = N_MBQUAD1(ksi, zeta);
223  base(gg, 2) = N_MBQUAD2(ksi, zeta);
224  base(gg, 3) = N_MBQUAD3(ksi, zeta);
225  diff_base(gg, 0) = diffN_MBQUAD0x(zeta);
226  diff_base(gg, 1) = diffN_MBQUAD0y(ksi);
227  diff_base(gg, 2) = diffN_MBQUAD1x(zeta);
228  diff_base(gg, 3) = diffN_MBQUAD1y(ksi);
229  diff_base(gg, 4) = diffN_MBQUAD2x(zeta);
230  diff_base(gg, 5) = diffN_MBQUAD2y(ksi);
231  diff_base(gg, 6) = diffN_MBQUAD3x(zeta);
232  diff_base(gg, 7) = diffN_MBQUAD3y(ksi);
233  }
235  };
236 
237  const auto type = numeredEntFiniteElementPtr->getEntType();
238 
239  if (rule >= 0) {
240  switch (type) {
241  case MBTRI:
242  CHKERR set_integration_pts_for_tri();
243  break;
244  case MBQUAD:
246  rule);
247  CHKERR calc_base_for_quad();
248  break;
249  default:
250  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
251  "Element type not implemented: %d", type);
252  }
253 
254  } else {
255  // If rule is negative, set user defined integration points
256  CHKERR setGaussPts(order_row, order_col, order_data);
257  const size_t nb_gauss_pts = gaussPts.size2();
258  auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
259  auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
260  if (nb_gauss_pts) {
261  switch (type) {
262  case MBTRI:
263  base.resize(nb_gauss_pts, 3, false);
264  diff_base.resize(3, 2, false);
265  CHKERR ShapeMBTRI(&*base.data().begin(), &gaussPts(0, 0),
266  &gaussPts(1, 0), nb_gauss_pts);
267  std::copy(
269  dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
270  break;
271  case MBQUAD:
272  CHKERR calc_base_for_quad();
273  break;
274  default:
275  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
276  "Element type not implemented: %d", type);
277  }
278  }
279  }
281 }
282 
286  // Get spaces order/base and sense of entities.
287 
289 
290  // H1
291  if (dataH1.spacesOnEntities[MBEDGE].test(H1)) {
292  CHKERR getEntitySense<MBEDGE>(dataH1);
293  CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
294  }
295  if (dataH1.spacesOnEntities[MBTRI].test(H1)) {
296  CHKERR getEntitySense<MBTRI>(dataH1);
297  CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
298  }
299  if (dataH1.spacesOnEntities[MBQUAD].test(H1)) {
300  CHKERR getEntitySense<MBQUAD>(dataH1);
301  CHKERR getEntityDataOrder<MBQUAD>(dataH1, H1);
302  }
303 
304  // Hcurl
305  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
306  CHKERR getEntitySense<MBEDGE>(dataHcurl);
307  CHKERR getEntityDataOrder<MBEDGE>(dataHcurl, HCURL);
308  dataHcurl.spacesOnEntities[MBEDGE].set(HCURL);
309  }
310  if (dataH1.spacesOnEntities[MBTRI].test(HCURL)) {
311  CHKERR getEntitySense<MBTRI>(dataHcurl);
312  CHKERR getEntityDataOrder<MBTRI>(dataHcurl, HCURL);
313  dataHcurl.spacesOnEntities[MBTRI].set(HCURL);
314  }
315 
316  // Hdiv
317  if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
318  CHKERR getEntitySense<MBTRI>(dataHdiv);
319  CHKERR getEntityDataOrder<MBTRI>(dataHdiv, HDIV);
320  dataHdiv.spacesOnEntities[MBTRI].set(HDIV);
321  }
322 
323  // L2
324  if (dataH1.spacesOnEntities[MBTRI].test(L2)) {
325  CHKERR getEntitySense<MBTRI>(dataL2);
326  CHKERR getEntityDataOrder<MBTRI>(dataL2, L2);
327  dataL2.spacesOnEntities[MBTRI].set(L2);
328  }
329 
331 }
332 
336 
337  const size_t nb_nodes =
338  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).size2();
339  double *shape_functions =
340  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
341  const size_t nb_gauss_pts = gaussPts.size2();
342  coordsAtGaussPts.resize(nb_gauss_pts, 3, false);
343  for (int gg = 0; gg != nb_gauss_pts; ++gg)
344  for (int dd = 0; dd != 3; ++dd)
346  nb_nodes, &shape_functions[nb_nodes * gg], 1, &coords[dd], 3);
347 
349 }
350 
353  // Check if field for high-order geometry is set and if it is set calculate
354  // higher-order normals and face tangent vectors.
355  if (dataPtr->get<FieldName_mi_tag>().find(meshPositionsFieldName) !=
356  dataPtr->get<FieldName_mi_tag>().end()) {
357 
358  const Field *field_struture =
360  BitFieldId id = field_struture->getId();
361 
362  if ((numeredEntFiniteElementPtr->getBitFieldIdData() & id).none())
363  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
364  "no MESH_NODE_POSITIONS in element data");
365 
366  // Calculate normal for high-order geometry
372 
373  } else if(numeredEntFiniteElementPtr->getEntType() == MBTRI) {
374  hoCoordsAtGaussPts.resize(0, 0, false);
375  normalsAtGaussPts.resize(0, 0, false);
376  tangentOneAtGaussPts.resize(0, 0, false);
377  tangentTwoAtGaussPts.resize(0, 0, false);
378  } else {
379  hoCoordsAtGaussPts.resize(coordsAtGaussPts.size1(),
380  coordsAtGaussPts.size2(), false);
382  }
383 
385 }
386 
387 } // namespace MoFEM
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:67
DataForcesAndSourcesCore & dataL2
int getMaxRowOrder() const
Get max order of approximation for field in rows.
MatrixDouble gaussPts
Matrix of integration points.
field with continuous normal traction
Definition: definitions.h:173
virtual moab::Interface & get_moab()=0
#define diffN_MBQUAD3x(y)
Definition: fem_tools.h:75
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
Provide data structure for (tensor) field approximation.The Field is intended to provide support for ...
DataForcesAndSourcesCore & dataH1
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
int npoints
Definition: quad.h:29
virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement()
Determine approximation space and order of base functions.
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_dcopy.c:11
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
virtual MoFEMErrorCode setIntegrationPts()
Set integration points.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
std::string meshPositionsFieldName
Name of the field with geometry.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define diffN_MBQUAD0x(y)
Definition: fem_tools.h:69
#define diffN_MBQUAD0y(x)
Definition: fem_tools.h:70
boost::shared_ptr< const FEDofEntity_multiIndex > dataPtr
Pointer to finite element data dofs.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
#define diffN_MBQUAD2y(x)
Definition: fem_tools.h:74
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
#define QUAD_2D_TABLE_SIZE
Definition: quad.h:174
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
MoFEMErrorCode getNodesFieldData(DataForcesAndSourcesCore &data, const std::string &field_name) const
Get data on nodes.
#define diffN_MBQUAD2x(y)
Definition: fem_tools.h:73
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:68
DataForcesAndSourcesCore & dataHcurl
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:65
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
const BitFieldId & getId() const
Get unique field id.
field with continuous tangents
Definition: definitions.h:172
static constexpr std::array< double, 8 > diffShapeFunMBQUADAtCenter
Definition: Tools.hpp:112
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:87
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:53
#define CHKERR
Inline error check.
Definition: definitions.h:596
virtual MoFEMErrorCode opRhs(DataForcesAndSourcesCore &data, const bool do_vertices, const bool do_edges, const bool do_quads, const bool do_tris, const bool do_tets, const bool do_prisms, const bool error_if_no_base=true)
MoFEMErrorCode getEntityFieldData(DataForcesAndSourcesCore &data, const std::string &field_name, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
virtual MoFEMErrorCode calculateCoordinatesAtGaussPts()
Calculate coordinate at integration points.
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:194
virtual MoFEMErrorCode calculateAreaAndNormalAtIntegrationPts()
Calculate element area and normal of the face at integration points.
MultiIndex Tag for field name.
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForQuad(MatrixDouble &pts, const int edge0, const int edge1)
Definition: Tools.cpp:486
DataForcesAndSourcesCore & dataHdiv
#define diffN_MBQUAD1x(y)
Definition: fem_tools.h:71
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:66
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
structure to get information form mofem into DataForcesAndSourcesCore
virtual MoFEMErrorCode calculateHoNormal()
Calculate normal on curved elements.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
continuous field
Definition: definitions.h:171
virtual const Field * get_field_structure(const std::string &name)=0
get field structure
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
#define diffN_MBQUAD3y(x)
Definition: fem_tools.h:76
virtual MoFEMErrorCode calculateAreaAndNormal()
Calculate element area and normal of the face.
MoFEMErrorCode getSpacesAndBaseOnEntities(DataForcesAndSourcesCore &data) const
Get field approximation space and base on entities.
int getMaxColOrder() const
Get max order of approximation for field in columns.
int getMaxDataOrder() const
Get max order of approximation for data fields.
#define diffN_MBQUAD1y(x)
Definition: fem_tools.h:72
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
field with C-1 continuity
Definition: definitions.h:174