v0.10.0
FlatPrismElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1 /** \file FlatPrismElementForcesAndSourcesCore.cpp
2 
3 \brief Implementation of flat prism 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(hoCoordsAtGaussPtsF3, nOrmals_at_GaussPtF3,
28  tAngent1_at_GaussPtF3, tAngent2_at_GaussPtF3,
29  hoCoordsAtGaussPtsF4, nOrmals_at_GaussPtF4,
30  tAngent1_at_GaussPtF4, tAngent2_at_GaussPtF4) {
32  boost::shared_ptr<BaseFunction>(new FlatPrismPolynomialBase());
33 }
34 
37 
38  if (numeredEntFiniteElementPtr->getEntType() != MBPRISM)
41 
45 
47  int num_nodes;
48  const EntityHandle *conn;
49  CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
50  {
51  coords.resize(num_nodes * 3, false);
52  CHKERR mField.get_moab().get_coords(conn, num_nodes,
53  &*coords.data().begin());
54 
55  double diff_n[6];
56  CHKERR ShapeDiffMBTRI(diff_n);
57  normal.resize(6, false);
58  CHKERR ShapeFaceNormalMBTRI(diff_n, &coords[0], &normal[0]);
59  CHKERR ShapeFaceNormalMBTRI(diff_n, &coords[9], &normal[3]);
60  aRea[0] = cblas_dnrm2(3, &normal[0], 1) * 0.5;
61  aRea[1] = cblas_dnrm2(3, &normal[3], 1) * 0.5;
62  }
63 
65 
66  // H1
67  if ((dataH1.spacesOnEntities[MBVERTEX]).test(H1)) {
68  CHKERR getEntitySense<MBEDGE>(dataH1);
69  CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
70  CHKERR getEntitySense<MBTRI>(dataH1);
71  CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
72  }
73 
74  // H1
75  if ((dataH1.spacesOnEntities[MBEDGE]).test(HCURL)) {
76  CHKERR getEntitySense<MBEDGE>(data_curl);
77  CHKERR getEntityDataOrder<MBEDGE>(data_curl, HCURL);
78  CHKERR getEntitySense<MBTRI>(data_curl);
79  CHKERR getEntityDataOrder<MBTRI>(data_curl, HCURL);
80  }
81 
82  // Hdiv
83  if ((dataH1.spacesOnEntities[MBTRI]).test(HDIV)) {
84  CHKERR getEntitySense<MBTRI>(data_div);
85  CHKERR getEntityDataOrder<MBTRI>(data_div, HDIV);
86  }
87 
88  // L2
89  if ((dataH1.spacesOnEntities[MBTRI]).test(L2)) {
90  CHKERR getEntitySense<MBTRI>(data_l2);
91  CHKERR getEntityDataOrder<MBTRI>(data_l2, L2);
92  }
93 
94  int order_data = getMaxDataOrder();
95  int order_row = getMaxRowOrder();
96  int order_col = getMaxColOrder();
97  int rule = getRule(order_row, order_col, order_data);
98  int nb_gauss_pts;
99  // int rule = getRule(order);
100  if (rule >= 0) {
101  if (rule < QUAD_2D_TABLE_SIZE) {
102  if (QUAD_2D_TABLE[rule]->dim != 2) {
103  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
104  }
105  if (QUAD_2D_TABLE[rule]->order < rule) {
106  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
107  "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
108  }
109  nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
110  gaussPts.resize(3, nb_gauss_pts, false);
111  cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
112  &gaussPts(0, 0), 1);
113  cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
114  &gaussPts(1, 0), 1);
115  cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
116  &gaussPts(2, 0), 1);
117  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
118  false);
119  double *shape_ptr =
120  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
121  cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
122  1);
123  } else {
124  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
125  "rule > quadrature order %d < %d", rule, QUAD_2D_TABLE_SIZE);
126  nb_gauss_pts = 0;
127  }
128  } else {
129  CHKERR setGaussPts(order_row, order_col, order_data);
130  nb_gauss_pts = gaussPts.size2();
131  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
132  false);
133  if (nb_gauss_pts) {
135  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin(),
136  &gaussPts(0, 0), &gaussPts(1, 0), nb_gauss_pts);
137  }
138  }
139  if (nb_gauss_pts == 0)
141 
142  {
143  coordsAtGaussPts.resize(nb_gauss_pts, 6, false);
144  for (int gg = 0; gg < nb_gauss_pts; gg++) {
145  for (int dd = 0; dd < 3; dd++) {
147  3, &dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, 0), 1,
148  &coords[dd], 3);
149  coordsAtGaussPts(gg, 3 + dd) = cblas_ddot(
150  3, &dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, 0), 1,
151  &coords[9 + dd], 3);
152  }
153  }
154  }
155 
156  for (int space = HCURL; space != LASTSPACE; ++space)
157  if (dataOnElement[space]) {
158  dataH1.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE) =
159  dataOnElement[H1]->dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE);
160  }
161 
162  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
163  if (dataH1.bAse.test(b)) {
164  switch (static_cast<FieldApproximationBase>(b)) {
167  if (dataH1.spacesOnEntities[MBVERTEX].test(H1)) {
168  CHKERR getElementPolynomialBase()->getValue(
169  gaussPts,
170  boost::shared_ptr<BaseFunctionCtx>(new FlatPrismPolynomialBaseCtx(
172  H1, static_cast<FieldApproximationBase>(b), NOBASE)));
173  }
174  if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
175  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
176  "Not yet implemented");
177  }
178  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
179  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
180  "Not yet implemented");
181  }
182  if (dataH1.spacesOnEntities[MBTET].test(L2)) {
183  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
184  "Not yet implemented");
185  }
186  break;
187  default:
188  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
189  "Not yet implemented");
190  }
191  }
192  }
193 
194  auto check_field = [&]() {
195  auto field_it =
197  if (field_it != fieldsPtr->get<FieldName_mi_tag>().end())
198  if ((numeredEntFiniteElementPtr->getBitFieldIdData() &
199  (*field_it)->getId())
200  .any())
201  return true;
202  return false;
203  };
204 
205  // Check if field meshPositionsFieldName exist
206  if (check_field()) {
207 
208  hoCoordsAtGaussPtsF3.resize(nb_gauss_pts, 3, false);
209  nOrmals_at_GaussPtF3.resize(nb_gauss_pts, 3, false);
210  tAngent1_at_GaussPtF3.resize(nb_gauss_pts, 3, false);
211  tAngent2_at_GaussPtF3.resize(nb_gauss_pts, 3, false);
212  hoCoordsAtGaussPtsF4.resize(nb_gauss_pts, 3, false);
213  nOrmals_at_GaussPtF4.resize(nb_gauss_pts, 3, false);
214  tAngent1_at_GaussPtF4.resize(nb_gauss_pts, 3, false);
215  tAngent2_at_GaussPtF4.resize(nb_gauss_pts, 3, false);
221  } else {
222  hoCoordsAtGaussPtsF3.resize(0, 0, false);
223  nOrmals_at_GaussPtF3.resize(0, 0, false);
224  tAngent1_at_GaussPtF3.resize(0, 0, false);
225  tAngent2_at_GaussPtF3.resize(0, 0, false);
226  hoCoordsAtGaussPtsF4.resize(0, 0, false);
227  nOrmals_at_GaussPtF4.resize(0, 0, false);
228  tAngent1_at_GaussPtF4.resize(0, 0, false);
229  tAngent2_at_GaussPtF4.resize(0, 0, false);
230  }
231 
232  if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
233  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented yet");
234  }
235 
236  // Iterate over operators
238 
240 }
241 
243  ForcesAndSourcesCore *ptr) {
245  if (!(ptrFE = dynamic_cast<FlatPrismElementForcesAndSourcesCore *>(ptr)))
246  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
247  "User operator and finite element do not work together");
249 }
250 
251 } // namespace MoFEM
int getMaxRowOrder() const
Get max order of approximation for field in rows.
Deprecated interface functions.
MatrixDouble gaussPts
Matrix of integration points.
field with continuous normal traction
Definition: definitions.h:179
virtual moab::Interface & get_moab()=0
const Field_multiIndex * fieldsPtr
raw pointer to fields container
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
virtual MoFEMErrorCode opRhs(DataForcesAndSourcesCore &data, const bool error_if_no_base=false)
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
MoFEMErrorCode createDataOnElement()
Create a entity data on element object.
DataForcesAndSourcesCore & dataH1
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
MoFEMErrorCode operator()()
function is run for every finite element
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
int npoints
Definition: quad.h:29
PetscErrorCode ShapeFaceNormalMBTRI(double *diffN, const double *coords, double *normal)
Definition: fem_tools.c:241
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_dcopy.c:11
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
FieldSpace in [ 0, LASTSPACE )
Definition: definitions.h:181
double cblas_dnrm2(const int N, const double *X, const int incX)
Definition: cblas_dnrm2.c:12
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
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
MoFEMErrorCode getNodesFieldData(DataForcesAndSourcesCore &data, const std::string &field_name) const
Get data on nodes.
Calculate base functions on tetrahedralFIXME: Need moab and mofem finite element structure to work (t...
const int dim
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:206
constexpr int order
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
Class used to pass element data to calculate base functions on flat prism.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
field with continuous tangents
Definition: definitions.h:178
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
#define CHKERR
Inline error check.
Definition: definitions.h:604
MoFEMErrorCode getEntityFieldData(DataForcesAndSourcesCore &data, const std::string &field_name, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:194
MultiIndex Tag for field name.
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
structure to get information form mofem into DataForcesAndSourcesCore
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
continuous field
Definition: definitions.h:177
MoFEMErrorCode getSpacesAndBaseOnEntities(DataForcesAndSourcesCore &data) const
Get field approximation space and base on entities.
const std::array< boost::shared_ptr< DataForcesAndSourcesCore >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
std::bitset< LASTBASE > bAse
bases on element
int getMaxColOrder() const
Get max order of approximation for field in columns.
int getMaxDataOrder() const
Get max order of approximation for data fields.
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
field with C-1 continuity
Definition: definitions.h:180