v0.13.1
Public Member Functions | Private Member Functions | Private Attributes | List of all members
MoFEM::TriPolynomialBase Struct Reference

Calculate base functions on triangle. More...

#include <src/approximation/TriPolynomialBase.hpp>

Inheritance diagram for MoFEM::TriPolynomialBase:
[legend]
Collaboration diagram for MoFEM::TriPolynomialBase:
[legend]

Public Member Functions

MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const
 
 TriPolynomialBase ()
 
 ~TriPolynomialBase ()
 
MoFEMErrorCode getValue (MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
 
- Public Member Functions inherited from MoFEM::BaseFunction
MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, MoFEM::UnknownInterface **iface) const
 
virtual MoFEMErrorCode getValue (MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
 
virtual MoFEMErrorCode getValue (MatrixDouble &pts_x, MatrixDouble &pts_t, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
 
- Public Member Functions inherited from MoFEM::BaseFunctionUnknownInterface
virtual MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
 
virtual ~BaseFunctionUnknownInterface ()=default
 
- Public Member Functions inherited from MoFEM::UnknownInterface
virtual MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
 
template<class IFACE >
MoFEMErrorCode registerInterface (bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()=default
 

Private Member Functions

MoFEMErrorCode getValueH1 (MatrixDouble &pts)
 
MoFEMErrorCode getValueH1AinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueH1BernsteinBezierBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueL2 (MatrixDouble &pts)
 
MoFEMErrorCode getValueL2AinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueL2BernsteinBezierBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHdiv (MatrixDouble &pts)
 
MoFEMErrorCode getValueHcurl (MatrixDouble &pts)
 
MoFEMErrorCode getValueHdivAinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHcurlAinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHdivDemkowiczBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHcurlDemkowiczBase (MatrixDouble &pts)
 

Private Attributes

EntPolynomialBaseCtxcTx
 
ublas::matrix< MatrixDoubleN_face_edge
 
ublas::vector< MatrixDoubleN_face_bubble
 
ublas::matrix< MatrixDoublediffN_face_edge
 
ublas::vector< MatrixDoublediffN_face_bubble
 

Additional Inherited Members

- Static Public Member Functions inherited from MoFEM::UnknownInterface
static MoFEMErrorCode getLibVersion (Version &version)
 Get library version. More...
 
static MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version)
 Get database major version. More...
 
static MoFEMErrorCode setFileVersion (moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
 Get database major version. More...
 
static MoFEMErrorCode getInterfaceVersion (Version &version)
 Get database major version. More...
 

Detailed Description

Calculate base functions on triangle.

Definition at line 18 of file TriPolynomialBase.hpp.

Constructor & Destructor Documentation

◆ TriPolynomialBase()

TriPolynomialBase::TriPolynomialBase ( )

Definition at line 10 of file TriPolynomialBase.cpp.

10{}

◆ ~TriPolynomialBase()

TriPolynomialBase::~TriPolynomialBase ( )

Definition at line 11 of file TriPolynomialBase.cpp.

11{}

Member Function Documentation

◆ getValue()

MoFEMErrorCode TriPolynomialBase::getValue ( MatrixDouble pts,
boost::shared_ptr< BaseFunctionCtx ctx_ptr 
)
virtual

Reimplemented from MoFEM::BaseFunction.

Definition at line 847 of file TriPolynomialBase.cpp.

848 {
850
852
853 int nb_gauss_pts = pts.size2();
854 if (!nb_gauss_pts) {
856 }
857
858 if (pts.size1() < 2)
859 SETERRQ(
860 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
861 "Wrong dimension of pts, should be at least 3 rows with coordinates");
862
863 const FieldApproximationBase base = cTx->bAse;
864 EntitiesFieldData &data = cTx->dAta;
865
867 if (cTx->copyNodeBase == LASTBASE) {
868 data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 3,
869 false);
871 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
872 &pts(0, 0), &pts(1, 0), nb_gauss_pts);
873 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(3, 2, false);
875 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
876 } else {
877 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
878 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
879 data.dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(base) =
880 data.dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(cTx->copyNodeBase);
881 }
882 if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
883 (unsigned int)nb_gauss_pts) {
884 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
885 "Base functions or nodes has wrong number of integration points "
886 "for base %s",
888 }
889 }
890 auto set_node_derivative_for_all_gauss_pts = [&]() {
892 // In linear geometry derivatives are constant,
893 // this in expense of efficiency makes implementation
894 // consistent between vertices and other types of entities
895 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(nb_gauss_pts, 6,
896 false);
897 for (int gg = 0; gg != nb_gauss_pts; ++gg)
898 std::copy(Tools::diffShapeFunMBTRI.begin(),
900 &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(gg, 0));
902 };
903
904 CHKERR set_node_derivative_for_all_gauss_pts();
905
906 switch (cTx->sPace) {
907 case H1:
908 CHKERR getValueH1(pts);
909 break;
910 case HDIV:
911 CHKERR getValueHdiv(pts);
912 break;
913 case HCURL:
915 break;
916 case L2:
917 CHKERR getValueL2(pts);
918 break;
919 default:
920 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
921 }
922
924}
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
static const char *const ApproximationBaseNames[]
Definition: definitions.h:72
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:194
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:182
Class used to pass element data to calculate base functions on tet,triangle,edge.
const FieldApproximationBase copyNodeBase
const FieldApproximationBase bAse
data structure for finite element entity
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:104
MoFEMErrorCode getValueH1(MatrixDouble &pts)
MoFEMErrorCode getValueL2(MatrixDouble &pts)
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
EntPolynomialBaseCtx * cTx
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ getValueH1()

MoFEMErrorCode TriPolynomialBase::getValueH1 ( MatrixDouble pts)
private

Definition at line 21 of file TriPolynomialBase.cpp.

21 {
23
24 switch (cTx->bAse) {
28 break;
31 break;
32 default:
33 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
34 }
35
37}
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)

◆ getValueH1AinsworthBase()

MoFEMErrorCode TriPolynomialBase::getValueH1AinsworthBase ( MatrixDouble pts)
private

Definition at line 253 of file TriPolynomialBase.cpp.

253 {
255
256 EntitiesFieldData &data = cTx->dAta;
257 const FieldApproximationBase base = cTx->bAse;
258
259 if (cTx->basePolynomialsType0 == NULL)
260 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
261 "Polynomial type not set");
262
263 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
264 double *diffL, const int dim) =
266
267 int nb_gauss_pts = pts.size2();
268
269 if (data.spacesOnEntities[MBEDGE].test(H1)) {
270 // edges
271 if (data.dataOnEntities[MBEDGE].size() != 3)
272 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
273 "expected size of data.dataOnEntities[MBEDGE] is 3");
274
275 int sense[3], order[3];
276 double *H1edgeN[3], *diffH1edgeN[3];
277 for (int ee = 0; ee < 3; ee++) {
278
279 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
280 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
281 "sense of the edge unknown");
282
283 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
284 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
285 int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getOrder());
286 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
287 false);
288 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
289 2 * nb_dofs, false);
290 H1edgeN[ee] = &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
291 diffH1edgeN[ee] =
292 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
293 }
295 sense, order,
296 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
297 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
298 H1edgeN, diffH1edgeN, nb_gauss_pts, base_polynomials);
299 }
300
301 if (data.spacesOnEntities[MBTRI].test(H1)) {
302 // face
303 if (data.dataOnEntities[MBTRI].size() != 1) {
304 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
305 "expected that size data.dataOnEntities[MBTRI] is one");
306 }
307
308 int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][0].getOrder());
309 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, nb_dofs,
310 false);
311 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
312 2 * nb_dofs, false);
313 const int face_nodes[] = {0, 1, 2};
315 face_nodes, data.dataOnEntities[MBTRI][0].getOrder(),
316 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
317 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
318 &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
319 &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
320 nb_gauss_pts, base_polynomials);
321 }
322
324}
static Index< 'p', 3 > p
const int dim
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
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:17
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:191
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types

◆ getValueH1BernsteinBezierBase()

MoFEMErrorCode TriPolynomialBase::getValueH1BernsteinBezierBase ( MatrixDouble pts)
private

Definition at line 40 of file TriPolynomialBase.cpp.

40 {
42 EntitiesFieldData &data = cTx->dAta;
43 const std::string &field_name = cTx->fieldName;
44 int nb_gauss_pts = pts.size2();
45
46 auto get_alpha = [field_name](auto &data) -> MatrixInt & {
47 auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
48 if (!ptr)
49 ptr.reset(new MatrixInt());
50 return *ptr;
51 };
52
53 auto get_base = [field_name](auto &data) -> MatrixDouble & {
54 auto &ptr = data.getBBNSharedPtr(field_name);
55 if (!ptr)
56 ptr.reset(new MatrixDouble());
57 return *ptr;
58 };
59
60 auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
61 auto &ptr = data.getBBDiffNSharedPtr(field_name);
62 if (!ptr)
63 ptr.reset(new MatrixDouble());
64 return *ptr;
65 };
66
67 auto get_alpha_by_name_ptr =
68 [](auto &data,
69 const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
70 return data.getBBAlphaIndicesSharedPtr(field_name);
71 };
72
73 auto get_base_by_name_ptr =
74 [](auto &data,
75 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
76 return data.getBBNSharedPtr(field_name);
77 };
78
79 auto get_diff_base_by_name_ptr =
80 [](auto &data,
81 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
82 return data.getBBDiffNSharedPtr(field_name);
83 };
84
85 auto get_alpha_by_order_ptr =
86 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
87 return data.getBBAlphaIndicesByOrderSharedPtr(o);
88 };
89
90 auto get_base_by_order_ptr =
91 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
92 return data.getBBNByOrderSharedPtr(o);
93 };
94
95 auto get_diff_base_by_order_ptr =
96 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
97 return data.getBBDiffNByOrderSharedPtr(o);
98 };
99
100 auto &vert_get_n = get_base(data.dataOnEntities[MBVERTEX][0]);
101 auto &vert_get_diff_n = get_diff_base(data.dataOnEntities[MBVERTEX][0]);
102 vert_get_n.resize(nb_gauss_pts, 3, false);
103 vert_get_diff_n.resize(nb_gauss_pts, 6, false);
104
105 if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
106 (unsigned int)nb_gauss_pts)
107 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
108 "Base functions or nodes has wrong number of integration points "
109 "for base %s",
111 auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
112
113 auto &vertex_alpha = get_alpha(data.dataOnEntities[MBVERTEX][0]);
114 vertex_alpha.resize(3, 3, false);
115 vertex_alpha.clear();
116 for (int n = 0; n != 3; ++n)
117 vertex_alpha(n, n) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n];
118
120 1, lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
121 &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &vert_get_n(0, 0),
122 &vert_get_diff_n(0, 0));
123 for (int n = 0; n != 3; ++n) {
124 const int f = boost::math::factorial<double>(
125 data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n]);
126 for (int g = 0; g != nb_gauss_pts; ++g) {
127 vert_get_n(g, n) *= f;
128 for (int d = 0; d != 2; ++d)
129 vert_get_diff_n(g, 2 * n + d) *= f;
130 }
131 }
132
133 // edges
134 if (data.spacesOnEntities[MBEDGE].test(H1)) {
135 if (data.dataOnEntities[MBEDGE].size() != 3)
136 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
137 "Wrong size of ent data");
138
139 constexpr int edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
140 for (int ee = 0; ee != 3; ++ee) {
141 auto &ent_data = data.dataOnEntities[MBEDGE][ee];
142
143 if (ent_data.getSense() == 0)
144 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
145 "Sense of the edge unknown");
146
147 const int sense = ent_data.getSense();
148 const int order = ent_data.getOrder();
149 const int nb_dofs = NBEDGE_H1(ent_data.getOrder());
150
151 if (nb_dofs) {
152 if (get_alpha_by_order_ptr(ent_data, order)) {
153 get_alpha_by_name_ptr(ent_data, field_name) =
154 get_alpha_by_order_ptr(ent_data, order);
155 get_base_by_name_ptr(ent_data, field_name) =
156 get_base_by_order_ptr(ent_data, order);
157 get_diff_base_by_name_ptr(ent_data, field_name) =
158 get_diff_base_by_order_ptr(ent_data, order);
159 } else {
160 auto &get_n = get_base(ent_data);
161 auto &get_diff_n = get_diff_base(ent_data);
162 get_n.resize(nb_gauss_pts, nb_dofs, false);
163 get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
164
165 auto &edge_alpha = get_alpha(data.dataOnEntities[MBEDGE][ee]);
166 edge_alpha.resize(nb_dofs, 3, false);
168 &edge_alpha(0, 0));
169 if (sense == -1) {
170 for (int i = 0; i != edge_alpha.size1(); ++i) {
171 int a = edge_alpha(i, edges_nodes[ee][0]);
172 edge_alpha(i, edges_nodes[ee][0]) =
173 edge_alpha(i, edges_nodes[ee][1]);
174 edge_alpha(i, edges_nodes[ee][1]) = a;
175 }
176 }
178 order, lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
179 &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &get_n(0, 0),
180 &get_diff_n(0, 0));
181
182 get_alpha_by_order_ptr(ent_data, order) =
183 get_alpha_by_name_ptr(ent_data, field_name);
184 get_base_by_order_ptr(ent_data, order) =
185 get_base_by_name_ptr(ent_data, field_name);
186 get_diff_base_by_order_ptr(ent_data, order) =
187 get_diff_base_by_name_ptr(ent_data, field_name);
188 }
189 }
190 }
191 } else {
192 for (int ee = 0; ee != 3; ++ee) {
193 auto &ent_data = data.dataOnEntities[MBEDGE][ee];
194 ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
195 auto &get_n = get_base(ent_data);
196 auto &get_diff_n = get_diff_base(ent_data);
197 get_n.resize(nb_gauss_pts, 0, false);
198 get_diff_n.resize(nb_gauss_pts, 0, false);
199 }
200 }
201
202 if (data.spacesOnEntities[MBTRI].test(H1)) {
203 if (data.dataOnEntities[MBTRI].size() != 1)
204 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
205 "Wrong size ent of ent data");
206
207 auto &ent_data = data.dataOnEntities[MBTRI][0];
208 const int order = ent_data.getOrder();
209 const int nb_dofs = NBFACETRI_H1(order);
210 if (get_alpha_by_order_ptr(ent_data, order)) {
211 get_alpha_by_name_ptr(ent_data, field_name) =
212 get_alpha_by_order_ptr(ent_data, order);
213 get_base_by_name_ptr(ent_data, field_name) =
214 get_base_by_order_ptr(ent_data, order);
215 get_diff_base_by_name_ptr(ent_data, field_name) =
216 get_diff_base_by_order_ptr(ent_data, order);
217 } else {
218
219 auto &get_n = get_base(ent_data);
220 auto &get_diff_n = get_diff_base(ent_data);
221 get_n.resize(nb_gauss_pts, nb_dofs, false);
222 get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
223 if (nb_dofs) {
224 auto &tri_alpha = get_alpha(ent_data);
225 tri_alpha.resize(nb_dofs, 3, false);
226
229 order, lambda.size1(), tri_alpha.size1(), &tri_alpha(0, 0),
230 &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &get_n(0, 0),
231 &get_diff_n(0, 0));
232
233 get_alpha_by_order_ptr(ent_data, order) =
234 get_alpha_by_name_ptr(ent_data, field_name);
235 get_base_by_order_ptr(ent_data, order) =
236 get_base_by_name_ptr(ent_data, field_name);
237 get_diff_base_by_order_ptr(ent_data, order) =
238 get_diff_base_by_name_ptr(ent_data, field_name);
239 }
240 }
241 } else {
242 auto &ent_data = data.dataOnEntities[MBTRI][0];
243 ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
244 auto &get_n = get_base(ent_data);
245 auto &get_diff_n = get_diff_base(ent_data);
246 get_n.resize(nb_gauss_pts, 0, false);
247 get_diff_n.resize(nb_gauss_pts, 0, false);
248 }
249
251}
static Index< 'o', 3 > o
constexpr double a
@ NOBASE
Definition: definitions.h:59
FTensor::Index< 'n', SPACE_DIM > n
constexpr double lambda
FTensor::Index< 'i', SPACE_DIM > i
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
auto f
Definition: HenckyOps.hpp:5
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasMatrix< int > MatrixInt
Definition: Types.hpp:76
constexpr auto field_name
constexpr double g
static MoFEMErrorCode baseFunctionsTri(const int N, const int gdim, const int n_alpha, const int *alpha, const double *lambda, const double *grad_lambda, double *base, double *grad_base)
static MoFEMErrorCode generateIndicesTriTri(const int N, int *alpha)
static MoFEMErrorCode generateIndicesEdgeTri(const int N[], int *alpha[])

◆ getValueHcurl()

MoFEMErrorCode TriPolynomialBase::getValueHcurl ( MatrixDouble pts)
private

Definition at line 828 of file TriPolynomialBase.cpp.

828 {
830
831 switch (cTx->bAse) {
835 break;
838 break;
839 default:
840 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
841 }
842
844}
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)

◆ getValueHcurlAinsworthBase()

MoFEMErrorCode TriPolynomialBase::getValueHcurlAinsworthBase ( MatrixDouble pts)
private

Definition at line 648 of file TriPolynomialBase.cpp.

648 {
650
651 EntitiesFieldData &data = cTx->dAta;
652 const FieldApproximationBase base = cTx->bAse;
653 if (data.dataOnEntities[MBTRI].size() != 1)
654 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
655
656 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
657 double *diffL, const int dim) =
659
660 int nb_gauss_pts = pts.size2();
661
662 // Calculation H-curl on triangle faces
663 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
664 if (data.dataOnEntities[MBEDGE].size() != 3) {
665 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
666 }
667 int sense[3], order[3];
668 double *hcurl_edge_n[3];
669 double *diff_hcurl_edge_n[3];
670 for (int ee = 0; ee < 3; ee++) {
671 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
672 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
673 "data inconsistency");
674 }
675 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
676 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
677 int nb_dofs = NBEDGE_AINSWORTH_HCURL(
678 data.dataOnEntities[MBEDGE][ee].getOrder());
679 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
680 3 * nb_dofs, false);
681 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(
682 nb_gauss_pts, 2 * 3 * nb_dofs, false);
683 hcurl_edge_n[ee] =
684 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
685 diff_hcurl_edge_n[ee] =
686 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
687 }
689 sense, order,
690 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
691 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
692 hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
693 } else {
694 for (int ee = 0; ee < 3; ee++) {
695 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
696 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
697 false);
698 }
699 }
700
701 if (data.spacesOnEntities[MBTRI].test(HCURL)) {
702
703 // cerr << data.dataOnEntities[MBVERTEX][0].getN(base) << endl;
704 // cerr << data.dataOnEntities[MBVERTEX][0].getDiffN(base) << endl;
705 //
706 // face
707 if (data.dataOnEntities[MBTRI].size() != 1) {
708 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
709 }
710 int order = data.dataOnEntities[MBTRI][0].getOrder();
711 int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order);
712 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
713 false);
714 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
715 3 * 2 * nb_dofs, false);
716 // cerr << data.dataOnEntities[MBVERTEX][0].getDiffN(base) << endl;
717 int face_nodes[] = {0, 1, 2};
719 face_nodes, order,
720 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
721 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
722 &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
723 &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
724 nb_gauss_pts, base_polynomials);
725 // cerr << data.dataOnEntities[MBTRI][0].getN(base) << endl;
726 } else {
727 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
728 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
729 }
730
732}
#define NBFACETRI_AINSWORTH_HCURL(P)
#define NBEDGE_AINSWORTH_HCURL(P)
MoFEMErrorCode Hcurl_Ainsworth_FaceFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face H-curl functions.
Definition: Hcurl.cpp:1249
MoFEMErrorCode Hcurl_Ainsworth_EdgeBaseFunctions_MBTET_ON_FACE(int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Edge based H-curl base functions on face.
Definition: Hcurl.cpp:237

◆ getValueHcurlDemkowiczBase()

MoFEMErrorCode TriPolynomialBase::getValueHcurlDemkowiczBase ( MatrixDouble pts)
private

Definition at line 735 of file TriPolynomialBase.cpp.

735 {
737
738 EntitiesFieldData &data = cTx->dAta;
739 const FieldApproximationBase base = cTx->bAse;
740
741 int nb_gauss_pts = pts.size2();
742
743 // Calculation H-curl on triangle faces
744 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
745
746 if (data.dataOnEntities[MBEDGE].size() != 3)
747 SETERRQ1(
748 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
749 "wrong number of data structures on edges, should be three but is %d",
750 data.dataOnEntities[MBEDGE].size());
751
752 int sense[3], order[3];
753 double *hcurl_edge_n[3];
754 double *diff_hcurl_edge_n[3];
755
756 for (int ee = 0; ee != 3; ++ee) {
757
758 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
759 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
760 "orientation (sense) of edge is not set");
761
762 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
763 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
764 int nb_dofs = NBEDGE_DEMKOWICZ_HCURL(
765 data.dataOnEntities[MBEDGE][ee].getOrder());
766 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
767 3 * nb_dofs, false);
768 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(
769 nb_gauss_pts, 2 * 3 * nb_dofs, false);
770
771 hcurl_edge_n[ee] =
772 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
773 diff_hcurl_edge_n[ee] =
774 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
775 }
776
778 sense, order,
779 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
780 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
781 hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
782
783 } else {
784
785 // No DOFs on faces, resize base function matrices, indicating that no
786 // dofs on them.
787 for (int ee = 0; ee != 3; ++ee) {
788 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
789 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
790 false);
791 }
792 }
793
794 if (data.spacesOnEntities[MBTRI].test(HCURL)) {
795
796 // face
797 if (data.dataOnEntities[MBTRI].size() != 1)
798 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
799 "No data struture to keep base functions on face");
800
801 int order = data.dataOnEntities[MBTRI][0].getOrder();
802 int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order);
803 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
804 false);
805 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
806 3 * 2 * nb_dofs, false);
807
808 int face_nodes[] = {0, 1, 2};
810 face_nodes, order,
811 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
812 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
813 &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
814 &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
815 nb_gauss_pts);
816
817 } else {
818
819 // No DOFs on faces, resize base function matrices, indicating that no
820 // dofs on them.
821 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
822 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
823 }
824
826}
#define NBEDGE_DEMKOWICZ_HCURL(P)
#define NBFACETRI_DEMKOWICZ_HCURL(P)
MoFEMErrorCode Hcurl_Demkowicz_EdgeBaseFunctions_MBTRI(int *sense, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Edge based H-curl base functions on teriangle.
Definition: Hcurl.cpp:2105
MoFEMErrorCode Hcurl_Demkowicz_FaceBaseFunctions_MBTRI(int *faces_nodes, int p, double *n, double *diff_n, double *phi, double *diff_phi, int nb_integration_pts)
Face base interior function.
Definition: Hcurl.cpp:2433

◆ getValueHdiv()

MoFEMErrorCode TriPolynomialBase::getValueHdiv ( MatrixDouble pts)
private

Definition at line 631 of file TriPolynomialBase.cpp.

631 {
633
634 switch (cTx->bAse) {
637 return getValueHdivAinsworthBase(pts);
639 return getValueHdivDemkowiczBase(pts);
640 default:
641 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
642 }
643
645}
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)

◆ getValueHdivAinsworthBase()

MoFEMErrorCode TriPolynomialBase::getValueHdivAinsworthBase ( MatrixDouble pts)
private

Definition at line 530 of file TriPolynomialBase.cpp.

530 {
532
533 EntitiesFieldData &data = cTx->dAta;
534 const FieldApproximationBase base = cTx->bAse;
535 if (cTx->basePolynomialsType0 == NULL)
536 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
537 "Polynomial type not set");
538 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
539 double *diffL, const int dim) =
541
542 int nb_gauss_pts = pts.size2();
543
544 double *PHI_f_e[3];
545 double *PHI_f;
546
547 N_face_edge.resize(1, 3, false);
548 N_face_bubble.resize(1, false);
549 int face_order = data.dataOnEntities[MBTRI][0].getOrder();
550 // three edges on face
551 for (int ee = 0; ee < 3; ee++) {
552 N_face_edge(0, ee).resize(
553 nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(face_order), false);
554 PHI_f_e[ee] = &((N_face_edge(0, ee))(0, 0));
555 }
556 N_face_bubble[0].resize(nb_gauss_pts,
557 3 * NBFACETRI_AINSWORTH_FACE_HDIV(face_order), false);
558 PHI_f = &*(N_face_bubble[0].data().begin());
559
560 int face_nodes[3] = {0, 1, 2};
562 face_nodes, face_order,
563 &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0), NULL, PHI_f_e, NULL,
564 nb_gauss_pts, 3, base_polynomials);
566 face_nodes, face_order,
567 &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0), NULL, PHI_f, NULL,
568 nb_gauss_pts, 3, base_polynomials);
569
570 // set shape functions into data structure
571 if (data.dataOnEntities[MBTRI].size() != 1) {
572 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
573 }
574 data.dataOnEntities[MBTRI][0].getN(base).resize(
575 nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_HDIV(face_order), false);
576 int col = 0;
577 for (int oo = 0; oo < face_order; oo++) {
578 for (int ee = 0; ee < 3; ee++) {
579 for (int dd = 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
580 dd < 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++, col++) {
581 for (int gg = 0; gg < nb_gauss_pts; gg++) {
582 data.dataOnEntities[MBTRI][0].getN(base)(gg, col) =
583 N_face_edge(0, ee)(gg, dd);
584 }
585 }
586 }
587 for (int dd = 3 * NBFACETRI_AINSWORTH_FACE_HDIV(oo);
588 dd < 3 * NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++, col++) {
589 for (int gg = 0; gg < nb_gauss_pts; gg++) {
590 data.dataOnEntities[MBTRI][0].getN(base)(gg, col) =
591 N_face_bubble[0](gg, dd);
592 }
593 }
594 }
595
597}
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
#define NBFACETRI_AINSWORTH_HDIV(P)
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
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f_e[3], double *diff_phi_f_e[3], int gdim, int nb, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth .
Definition: Hdiv.cpp:37
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face bubble functions by Ainsworth .
Definition: Hdiv.cpp:160
ublas::matrix< MatrixDouble > N_face_edge
ublas::vector< MatrixDouble > N_face_bubble

◆ getValueHdivDemkowiczBase()

MoFEMErrorCode TriPolynomialBase::getValueHdivDemkowiczBase ( MatrixDouble pts)
private

Definition at line 599 of file TriPolynomialBase.cpp.

599 {
601
602 EntitiesFieldData &data = cTx->dAta;
603 // set shape functions into data structure
604 if (data.dataOnEntities[MBTRI].size() != 1) {
605 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
606 }
607 const FieldApproximationBase base = cTx->bAse;
608 if (base != DEMKOWICZ_JACOBI_BASE) {
609 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
610 "This should be used only with DEMKOWICZ_JACOBI_BASE "
611 "but base is %s",
613 }
614 int order = data.dataOnEntities[MBTRI][0].getOrder();
615 int nb_gauss_pts = pts.size2();
616 data.dataOnEntities[MBTRI][0].getN(base).resize(
617 nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
618 double *phi_f = &*data.dataOnEntities[MBTRI][0].getN(base).data().begin();
621 int face_nodes[3] = {0, 1, 2};
623 face_nodes, order,
624 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
625 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(), phi_f,
626 NULL, nb_gauss_pts, 3);
627
629}
#define NBFACETRI_DEMKOWICZ_HDIV(P)
MoFEMErrorCode Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb)
Definition: Hdiv.cpp:617

◆ getValueL2()

MoFEMErrorCode TriPolynomialBase::getValueL2 ( MatrixDouble pts)
private

Definition at line 326 of file TriPolynomialBase.cpp.

326 {
328
329 switch (cTx->bAse) {
333 break;
336 break;
337 default:
338 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
339 }
340
342}
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)

◆ getValueL2AinsworthBase()

MoFEMErrorCode TriPolynomialBase::getValueL2AinsworthBase ( MatrixDouble pts)
private

Definition at line 344 of file TriPolynomialBase.cpp.

344 {
346
347 EntitiesFieldData &data = cTx->dAta;
348 const FieldApproximationBase base = cTx->bAse;
349 if (cTx->basePolynomialsType0 == NULL)
350 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
351 "Polynomial type not set");
352 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
353 double *diffL, const int dim) =
355
356 int nb_gauss_pts = pts.size2();
357
358 data.dataOnEntities[MBTRI][0].getN(base).resize(
359 nb_gauss_pts, NBFACETRI_L2(data.dataOnEntities[MBTRI][0].getOrder()),
360 false);
361 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(
362 nb_gauss_pts,
363 2 * NBFACETRI_L2(data.dataOnEntities[MBTRI][0].getOrder()), false);
364
366 data.dataOnEntities[MBTRI][0].getOrder(),
367 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
368 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
369 &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
370 &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
371 nb_gauss_pts, base_polynomials);
372
374}
#define NBFACETRI_L2(P)
Number of base functions on triangle for L2 space.
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTRI(int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Get base functions on triangle for L2 space.
Definition: l2.c:19

◆ getValueL2BernsteinBezierBase()

MoFEMErrorCode TriPolynomialBase::getValueL2BernsteinBezierBase ( MatrixDouble pts)
private

Definition at line 377 of file TriPolynomialBase.cpp.

377 {
379
380 EntitiesFieldData &data = cTx->dAta;
381 const std::string &field_name = cTx->fieldName;
382 int nb_gauss_pts = pts.size2();
383
384 auto get_alpha = [field_name](auto &data) -> MatrixInt & {
385 auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
386 if (!ptr)
387 ptr.reset(new MatrixInt());
388 return *ptr;
389 };
390
391 auto get_base = [field_name](auto &data) -> MatrixDouble & {
392 auto &ptr = data.getBBNSharedPtr(field_name);
393 if (!ptr)
394 ptr.reset(new MatrixDouble());
395 return *ptr;
396 };
397
398 auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
399 auto &ptr = data.getBBDiffNSharedPtr(field_name);
400 if (!ptr)
401 ptr.reset(new MatrixDouble());
402 return *ptr;
403 };
404
405 auto get_alpha_by_name_ptr =
406 [](auto &data,
407 const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
408 return data.getBBAlphaIndicesSharedPtr(field_name);
409 };
410
411 auto get_base_by_name_ptr =
412 [](auto &data,
413 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
414 return data.getBBNSharedPtr(field_name);
415 };
416
417 auto get_diff_base_by_name_ptr =
418 [](auto &data,
419 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
420 return data.getBBDiffNSharedPtr(field_name);
421 };
422
423 auto get_alpha_by_order_ptr =
424 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
425 return data.getBBAlphaIndicesByOrderSharedPtr(o);
426 };
427
428 auto get_base_by_order_ptr =
429 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
430 return data.getBBNByOrderSharedPtr(o);
431 };
432
433 auto get_diff_base_by_order_ptr =
434 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
435 return data.getBBDiffNByOrderSharedPtr(o);
436 };
437
438 if (data.spacesOnEntities[MBTRI].test(L2)) {
439 if (data.dataOnEntities[MBTRI].size() != 1)
440 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
441 "Wrong size ent of ent data");
442
443 if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
444 (unsigned int)nb_gauss_pts)
445 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
446 "Base functions or nodes has wrong number of integration points "
447 "for base %s",
449 auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
450
451 auto &ent_data = data.dataOnEntities[MBTRI][0];
452 const int order = ent_data.getOrder();
453 const int nb_dofs = NBFACETRI_L2(order);
454
455 if (get_alpha_by_order_ptr(ent_data, order)) {
456 get_alpha_by_name_ptr(ent_data, field_name) =
457 get_alpha_by_order_ptr(ent_data, order);
458 get_base_by_name_ptr(ent_data, field_name) =
459 get_base_by_order_ptr(ent_data, order);
460 get_diff_base_by_name_ptr(ent_data, field_name) =
461 get_diff_base_by_order_ptr(ent_data, order);
462 } else {
463
464 auto &get_n = get_base(ent_data);
465 auto &get_diff_n = get_diff_base(ent_data);
466 get_n.resize(nb_gauss_pts, nb_dofs, false);
467 get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
468 if (nb_dofs) {
469
470 if (order == 0) {
471
472 if (nb_dofs != 1)
473 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
474 "Inconsistent number of DOFs");
475
476 auto &tri_alpha = get_alpha(ent_data);
477 tri_alpha.clear();
478 get_n(0, 0) = 1;
479 get_diff_n.clear();
480
481 } else {
482
483 if (nb_dofs != 3 + 3 * NBEDGE_H1(order) + NBFACETRI_H1(order))
484 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
485 "Inconsistent number of DOFs");
486
487 auto &tri_alpha = get_alpha(ent_data);
488 tri_alpha.resize(nb_dofs, 3, false);
489
491 &tri_alpha(0, 0));
492
493 if (order > 1) {
494 std::array<int, 3> face_n{order, order, order};
495 std::array<int *, 3> face_edge_ptr{
496 &tri_alpha(3, 0), &tri_alpha(3 + NBEDGE_H1(order), 0),
497 &tri_alpha(3 + 2 * NBEDGE_H1(order), 0)};
499 face_n.data(), face_edge_ptr.data());
500 if (order > 2)
502 order, &tri_alpha(3 + 3 * NBEDGE_H1(order), 0));
503 }
505 order, lambda.size1(), tri_alpha.size1(), &tri_alpha(0, 0),
506 &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &get_n(0, 0),
507 &get_diff_n(0, 0));
508 }
509
510 get_alpha_by_order_ptr(ent_data, order) =
511 get_alpha_by_name_ptr(ent_data, field_name);
512 get_base_by_order_ptr(ent_data, order) =
513 get_base_by_name_ptr(ent_data, field_name);
514 get_diff_base_by_order_ptr(ent_data, order) =
515 get_diff_base_by_name_ptr(ent_data, field_name);
516 }
517 }
518 } else {
519 auto &ent_data = data.dataOnEntities[MBTRI][0];
520 ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
521 auto &get_n = get_base(ent_data);
522 auto &get_diff_n = get_diff_base(ent_data);
523 get_n.resize(nb_gauss_pts, 0, false);
524 get_diff_n.resize(nb_gauss_pts, 0, false);
525 }
526
528}
static MoFEMErrorCode generateIndicesVertexTri(const int N, int *alpha)

◆ query_interface()

MoFEMErrorCode TriPolynomialBase::query_interface ( boost::typeindex::type_index  type_index,
UnknownInterface **  iface 
) const
virtual

Reimplemented from MoFEM::BaseFunction.

Definition at line 14 of file TriPolynomialBase.cpp.

15 {
17 *iface = const_cast<TriPolynomialBase *>(this);
19}
Calculate base functions on triangle.

Member Data Documentation

◆ cTx

EntPolynomialBaseCtx* MoFEM::TriPolynomialBase::cTx
private

Definition at line 29 of file TriPolynomialBase.hpp.

◆ diffN_face_bubble

ublas::vector<MatrixDouble> MoFEM::TriPolynomialBase::diffN_face_bubble
private

Definition at line 42 of file TriPolynomialBase.hpp.

◆ diffN_face_edge

ublas::matrix<MatrixDouble> MoFEM::TriPolynomialBase::diffN_face_edge
private

Definition at line 41 of file TriPolynomialBase.hpp.

◆ N_face_bubble

ublas::vector<MatrixDouble> MoFEM::TriPolynomialBase::N_face_bubble
private

Definition at line 40 of file TriPolynomialBase.hpp.

◆ N_face_edge

ublas::matrix<MatrixDouble> MoFEM::TriPolynomialBase::N_face_edge
private

Definition at line 39 of file TriPolynomialBase.hpp.


The documentation for this struct was generated from the following files: