v0.9.0
FlatPrismPolynomialBase.cpp
Go to the documentation of this file.
1 /** \file FlatPrismPolynomialBase.cpp
2 \brief Implementation of Ainsworth-Cole H1 base on edge
3 */
4 
5 /* This file is part of MoFEM.
6 * MoFEM is free software: you can redistribute it and/or modify it under
7 * the terms of the GNU Lesser General Public License as published by the
8 * Free Software Foundation, either version 3 of the License, or (at your
9 * option) any later version.
10 *
11 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14 * License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
18 
19 
20 using namespace MoFEM;
21 
23  const MOFEMuuid &uuid, BaseFunctionUnknownInterface **iface) const {
24 
26  *iface = NULL;
27  if (uuid == IDD_FLATPRISM_BASE_FUNCTION) {
28  *iface = const_cast<FlatPrismPolynomialBaseCtx *>(this);
30  } else {
31  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "wrong interference");
32  }
35 }
36 
38  DataForcesAndSourcesCore &data, moab::Interface &moab,
39  const NumeredEntFiniteElement *fe_ptr, const FieldSpace space,
40  const FieldApproximationBase base,
41  const FieldApproximationBase copy_node_base)
42  : EntPolynomialBaseCtx(data, space, base, copy_node_base), mOab(moab),
43  fePtr(fe_ptr) {
44  CHKERR setBase();
45  CHKERRABORT(PETSC_COMM_WORLD, ierr);
46 }
48 
50  const MOFEMuuid &uuid, BaseFunctionUnknownInterface **iface) const {
51 
53  *iface = NULL;
54  if (uuid == IDD_FLATPRISM_BASE_FUNCTION) {
55  *iface = const_cast<FlatPrismPolynomialBase *>(this);
57  } else {
58  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "wrong interference");
59  }
62 }
63 
66 
69  boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
71 
73  CHKERR ctx_ptr->query_interface(IDD_FLATPRISM_BASE_FUNCTION, &iface);
74  cTx = reinterpret_cast<FlatPrismPolynomialBaseCtx *>(iface);
75  if (!cTx->fePtr)
76  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
77  "Pointer to element should be given "
78  "when EntPolynomialBaseCtx is constructed "
79  "(use different constructor)");
80 
81  int nb_gauss_pts = pts.size2();
82  if (!nb_gauss_pts)
84 
85  if (pts.size1() < 1)
86  SETERRQ(
87  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
88  "Wrong dimension of pts, should be at least 3 rows with coordinates");
89 
90  const FieldApproximationBase base = cTx->bAse;
92 
93  if (cTx->copyNodeBase == LASTBASE)
94  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
95  else
96  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
97  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
98 
99  data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 6, false);
100  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(nb_gauss_pts, 12,
101  false);
102  if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
103  (unsigned int)nb_gauss_pts) {
104  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
105  "Base functions or nodes has wrong number of integration points "
106  "for base %s",
107  ApproximationBaseNames[base]);
108  }
109 
110  data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 6, false);
111  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(nb_gauss_pts, 12,
112  false);
113  N.resize(nb_gauss_pts, 3, false);
114  diffN.resize(3, 2, false);
115  CHKERR ShapeMBTRI(&*N.data().begin(), &pts(0, 0), &pts(1, 0), nb_gauss_pts);
116  std::copy(Tools::diffShapeFunMBTRI.begin(), Tools::diffShapeFunMBTRI.end(),
117  &*diffN.data().begin());
118 
119  // This is needed to have proper order of nodes on faces
120  CHKERR cTx->mOab.get_connectivity(cTx->fePtr->getEnt(), connPrism, numNodes,
121  true);
122  SideNumber_multiIndex &side_table =
123  const_cast<SideNumber_multiIndex &>(cTx->fePtr->getSideNumberTable());
124  SideNumber_multiIndex::nth_index<1>::type::iterator siit3 =
125  side_table.get<1>().find(boost::make_tuple(MBTRI, 3));
126  SideNumber_multiIndex::nth_index<1>::type::iterator siit4 =
127  side_table.get<1>().find(boost::make_tuple(MBTRI, 4));
128  if (siit3 == side_table.get<1>().end())
129  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
130  if (siit4 == side_table.get<1>().end())
131  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
132  CHKERR cTx->mOab.get_connectivity(siit3->get()->ent, connFace3, numNodes,
133  true);
134  CHKERR cTx->mOab.get_connectivity(siit4->get()->ent, connFace4, numNodes,
135  true);
136 
137  for (int nn = 0; nn < 3; nn++) {
138  faceNodes[0][nn] = std::distance(
139  connPrism, std::find(connPrism, connPrism + 3, connFace3[nn]));
140  faceNodes[1][nn] = std::distance(
141  connPrism + 3, std::find(connPrism + 3, connPrism + 6, connFace4[nn]));
142  for (int gg = 0; gg < nb_gauss_pts; gg++) {
143  double val = N(gg, nn);
144  double val_x = diffN(nn, 0);
145  double val_y = diffN(nn, 1);
146  data.dataOnEntities[MBVERTEX][0].getN(base)(gg, nn) = val;
147  data.dataOnEntities[MBVERTEX][0].getDiffN(base)(gg, 2 * nn + 0) = val_x;
148  data.dataOnEntities[MBVERTEX][0].getDiffN(base)(gg, 2 * nn + 1) = val_y;
149  data.dataOnEntities[MBVERTEX][0].getN(base)(gg, 3 + nn) = val;
150  data.dataOnEntities[MBVERTEX][0].getDiffN(base)(gg, 6 + 2 * nn + 0) =
151  val_x;
152  data.dataOnEntities[MBVERTEX][0].getDiffN(base)(gg, 6 + 2 * nn + 1) =
153  val_y;
154  }
155  }
156 
157  switch (cTx->sPace) {
158  case H1:
159  CHKERR getValueH1(pts);
160  break;
161  case HDIV:
162  CHKERR getValueHdiv(pts);
163  break;
164  case HCURL:
165  CHKERR getValueHcurl(pts);
166  break;
167  case L2:
168  CHKERR getValueL2(pts);
169  break;
170  default:
171  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
172  }
173 
175 }
176 
179 
181  const FieldApproximationBase base = cTx->bAse;
182  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
183  double *diffL, const int dim) =
185 
186  int nb_gauss_pts = pts.size2();
187 
188  // edges
189  if (data.dataOnEntities[MBEDGE].size() != 9)
190  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
191 
192  int sense[9], order[9];
193  double *H1edgeN[9], *diffH1edgeN[9];
194 
195  auto set_edge_base_data = [&](const int ee) {
197  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
198  if (ent_data.getSense() == 0)
199  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
200  sense[ee] = ent_data.getSense();
201  order[ee] = ent_data.getDataOrder();
202  int nb_dofs = NBEDGE_H1(ent_data.getDataOrder());
203  ent_data.getN(base).resize(nb_gauss_pts, nb_dofs, false);
204  ent_data.getDiffN(base).resize(nb_gauss_pts, 2 * nb_dofs, false);
205  H1edgeN[ee] = &*ent_data.getN(base).data().begin();
206  diffH1edgeN[ee] = &*ent_data.getDiffN(base).data().begin();
208  };
209 
210  if ((data.spacesOnEntities[MBEDGE]).test(H1)) {
211 
212  for (int ee = 0; ee != 3; ++ee)
213  CHKERR set_edge_base_data(ee);
214  for (int ee = 6; ee != 9; ++ee)
215  CHKERR set_edge_base_data(ee);
216 
217  // shape functions on face 3
219  &sense[0], &order[0], &*N.data().begin(), &*diffN.data().begin(),
220  &H1edgeN[0], &diffH1edgeN[0], nb_gauss_pts, base_polynomials);
221  // shape functions on face 4
223  &sense[6], &order[6], &*N.data().begin(), &*diffN.data().begin(),
224  &H1edgeN[6], &diffH1edgeN[6], nb_gauss_pts, base_polynomials);
225  }
226 
227  // face
228  if (data.dataOnEntities[MBTRI].size() != 5)
229  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
230 
231  if ((data.spacesOnEntities[MBTRI]).test(H1)) {
232  for (int ff = 3; ff <= 4; ff++) {
233  int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][ff].getDataOrder());
234  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
235  false);
236  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
237  2 * nb_dofs, false);
239  faceNodes[ff - 3], data.dataOnEntities[MBTRI][ff].getDataOrder(),
240  &*N.data().begin(), &*diffN.data().begin(),
241  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin(),
242  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin(),
243  nb_gauss_pts, base_polynomials);
244  }
245  }
246 
248 }
249 
252  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
254 }
255 
257  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
258 }
259 
261  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
262 }
MoFEM interface unique ID.
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, BaseFunctionUnknownInterface **iface) const
field with continuous normal traction
Definition: definitions.h:173
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
static const MOFEMuuid IDD_FLATPRISM_BASE_FUNCTION
PetscErrorCode H1_EdgeShapeFunctions_MBTRI(int *sense, int *p, double *N, double *diffN, double *edgeN[3], double *diff_edgeN[3], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
H1_EdgeShapeFunctions_MBTRI.
Definition: h1.c:27
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
Partitioned (Indexed) Finite Element in Problem.
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, MoFEM::BaseFunctionUnknownInterface **iface) const
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
Class used to pass element data to calculate base functions on tet,triangle,edge.
const FieldApproximationBase bAse
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, BaseFunctionUnknownInterface **iface) const
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
PetscErrorCode H1_FaceShapeFunctions_MBTRI(const int *face_nodes, int p, double *N, double *diffN, double *faceN, double *diff_faceN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:162
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
const FieldApproximationBase copyNodeBase
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
static const char *const ApproximationBaseNames[]
Definition: definitions.h:158
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
FieldApproximationBase
approximation base
Definition: definitions.h:144
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, BaseFunctionUnknownInterface **iface) const
DataForcesAndSourcesCore & dAta
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
field with continuous tangents
Definition: definitions.h:172
const NumeredEntFiniteElement * fePtr
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
FieldSpace
approximation spaces
Definition: definitions.h:168
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:87
SideNumber_multiIndex & getSideNumberTable() const
#define CHKERR
Inline error check.
Definition: definitions.h:596
FlatPrismPolynomialBaseCtx * cTx
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:194
FlatPrismPolynomialBaseCtx(DataForcesAndSourcesCore &data, moab::Interface &moab, const NumeredEntFiniteElement *fe_ptr, const FieldSpace space, const FieldApproximationBase base, const FieldApproximationBase copy_node_base=LASTBASE)
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
#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
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< hashed_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, char, &SideNumber::side_number > > >, ordered_non_unique< const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
EntityHandle getEnt() const
Get element entity.
MoFEMErrorCode getValueL2(MatrixDouble &pts)
MoFEMErrorCode getValueH1(MatrixDouble &pts)
field with C-1 continuity
Definition: definitions.h:174