v0.13.2
Loading...
Searching...
No Matches
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 ()=default
 
virtual ~TriPolynomialBase ()=default
 
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 16 of file TriPolynomialBase.hpp.

Constructor & Destructor Documentation

◆ TriPolynomialBase()

MoFEM::TriPolynomialBase::TriPolynomialBase ( )
default

◆ ~TriPolynomialBase()

virtual MoFEM::TriPolynomialBase::~TriPolynomialBase ( )
virtualdefault

Member Function Documentation

◆ getValue()

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

Reimplemented from MoFEM::BaseFunction.

Definition at line 842 of file TriPolynomialBase.cpp.

843 {
845
847
848 int nb_gauss_pts = pts.size2();
849 if (!nb_gauss_pts) {
851 }
852
853 if (pts.size1() < 2)
854 SETERRQ(
855 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
856 "Wrong dimension of pts, should be at least 3 rows with coordinates");
857
858 const FieldApproximationBase base = cTx->bAse;
859 EntitiesFieldData &data = cTx->dAta;
860
862 if (cTx->copyNodeBase == LASTBASE) {
863 data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 3,
864 false);
866 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
867 &pts(0, 0), &pts(1, 0), nb_gauss_pts);
868 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(3, 2, false);
870 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
871 } else {
872 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
873 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
874 data.dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(base) =
875 data.dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(cTx->copyNodeBase);
876 }
877 if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
878 (unsigned int)nb_gauss_pts) {
879 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
880 "Base functions or nodes has wrong number of integration points "
881 "for base %s",
883 }
884 }
885 auto set_node_derivative_for_all_gauss_pts = [&]() {
887 // In linear geometry derivatives are constant,
888 // this in expense of efficiency makes implementation
889 // consistent between vertices and other types of entities
890 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(nb_gauss_pts, 6,
891 false);
892 for (int gg = 0; gg != nb_gauss_pts; ++gg)
893 std::copy(Tools::diffShapeFunMBTRI.begin(),
895 &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(gg, 0));
897 };
898
899 CHKERR set_node_derivative_for_all_gauss_pts();
900
901 switch (cTx->sPace) {
902 case H1:
903 CHKERR getValueH1(pts);
904 break;
905 case HDIV:
906 CHKERR getValueHdiv(pts);
907 break;
908 case HCURL:
910 break;
911 case L2:
912 CHKERR getValueL2(pts);
913 break;
914 default:
915 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
916 }
917
919}
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 16 of file TriPolynomialBase.cpp.

16 {
18
19 switch (cTx->bAse) {
23 break;
26 break;
27 default:
28 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
29 }
30
32}
@ 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 248 of file TriPolynomialBase.cpp.

248 {
250
251 EntitiesFieldData &data = cTx->dAta;
252 const FieldApproximationBase base = cTx->bAse;
253
254 if (cTx->basePolynomialsType0 == NULL)
255 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
256 "Polynomial type not set");
257
258 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
259 double *diffL, const int dim) =
261
262 int nb_gauss_pts = pts.size2();
263
264 if (data.spacesOnEntities[MBEDGE].test(H1)) {
265 // edges
266 if (data.dataOnEntities[MBEDGE].size() != 3)
267 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
268 "expected size of data.dataOnEntities[MBEDGE] is 3");
269
270 int sense[3], order[3];
271 double *H1edgeN[3], *diffH1edgeN[3];
272 for (int ee = 0; ee < 3; ee++) {
273
274 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
275 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
276 "sense of the edge unknown");
277
278 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
279 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
280 int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getOrder());
281 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
282 false);
283 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
284 2 * nb_dofs, false);
285 H1edgeN[ee] = &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
286 diffH1edgeN[ee] =
287 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
288 }
290 sense, order,
291 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
292 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
293 H1edgeN, diffH1edgeN, nb_gauss_pts, base_polynomials);
294 }
295
296 if (data.spacesOnEntities[MBTRI].test(H1)) {
297 // face
298 if (data.dataOnEntities[MBTRI].size() != 1) {
299 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
300 "expected that size data.dataOnEntities[MBTRI] is one");
301 }
302
303 int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][0].getOrder());
304 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, nb_dofs,
305 false);
306 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
307 2 * nb_dofs, false);
308 const int face_nodes[] = {0, 1, 2};
310 face_nodes, data.dataOnEntities[MBTRI][0].getOrder(),
311 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
312 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
313 &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
314 &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
315 nb_gauss_pts, base_polynomials);
316 }
317
319}
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 35 of file TriPolynomialBase.cpp.

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

823 {
825
826 switch (cTx->bAse) {
830 break;
833 break;
834 default:
835 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
836 }
837
839}
@ 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 643 of file TriPolynomialBase.cpp.

643 {
645
646 EntitiesFieldData &data = cTx->dAta;
647 const FieldApproximationBase base = cTx->bAse;
648 if (data.dataOnEntities[MBTRI].size() != 1)
649 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
650
651 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
652 double *diffL, const int dim) =
654
655 int nb_gauss_pts = pts.size2();
656
657 // Calculation H-curl on triangle faces
658 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
659 if (data.dataOnEntities[MBEDGE].size() != 3) {
660 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
661 }
662 int sense[3], order[3];
663 double *hcurl_edge_n[3];
664 double *diff_hcurl_edge_n[3];
665 for (int ee = 0; ee < 3; ee++) {
666 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
667 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
668 "data inconsistency");
669 }
670 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
671 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
672 int nb_dofs = NBEDGE_AINSWORTH_HCURL(
673 data.dataOnEntities[MBEDGE][ee].getOrder());
674 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
675 3 * nb_dofs, false);
676 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(
677 nb_gauss_pts, 2 * 3 * nb_dofs, false);
678 hcurl_edge_n[ee] =
679 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
680 diff_hcurl_edge_n[ee] =
681 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
682 }
684 sense, order,
685 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
686 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
687 hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
688 } else {
689 for (int ee = 0; ee < 3; ee++) {
690 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
691 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
692 false);
693 }
694 }
695
696 if (data.spacesOnEntities[MBTRI].test(HCURL)) {
697
698 // cerr << data.dataOnEntities[MBVERTEX][0].getN(base) << endl;
699 // cerr << data.dataOnEntities[MBVERTEX][0].getDiffN(base) << endl;
700 //
701 // face
702 if (data.dataOnEntities[MBTRI].size() != 1) {
703 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
704 }
705 int order = data.dataOnEntities[MBTRI][0].getOrder();
706 int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order);
707 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
708 false);
709 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
710 3 * 2 * nb_dofs, false);
711 // cerr << data.dataOnEntities[MBVERTEX][0].getDiffN(base) << endl;
712 int face_nodes[] = {0, 1, 2};
714 face_nodes, order,
715 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
716 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
717 &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
718 &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
719 nb_gauss_pts, base_polynomials);
720 // cerr << data.dataOnEntities[MBTRI][0].getN(base) << endl;
721 } else {
722 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
723 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
724 }
725
727}
#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 730 of file TriPolynomialBase.cpp.

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

626 {
628
629 switch (cTx->bAse) {
632 return getValueHdivAinsworthBase(pts);
634 return getValueHdivDemkowiczBase(pts);
635 default:
636 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
637 }
638
640}
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)

◆ getValueHdivAinsworthBase()

MoFEMErrorCode TriPolynomialBase::getValueHdivAinsworthBase ( MatrixDouble pts)
private

Definition at line 525 of file TriPolynomialBase.cpp.

525 {
527
528 EntitiesFieldData &data = cTx->dAta;
529 const FieldApproximationBase base = cTx->bAse;
530 if (cTx->basePolynomialsType0 == NULL)
531 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
532 "Polynomial type not set");
533 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
534 double *diffL, const int dim) =
536
537 int nb_gauss_pts = pts.size2();
538
539 double *PHI_f_e[3];
540 double *PHI_f;
541
542 N_face_edge.resize(1, 3, false);
543 N_face_bubble.resize(1, false);
544 int face_order = data.dataOnEntities[MBTRI][0].getOrder();
545 // three edges on face
546 for (int ee = 0; ee < 3; ee++) {
547 N_face_edge(0, ee).resize(
548 nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(face_order), false);
549 PHI_f_e[ee] = &((N_face_edge(0, ee))(0, 0));
550 }
551 N_face_bubble[0].resize(nb_gauss_pts,
552 3 * NBFACETRI_AINSWORTH_FACE_HDIV(face_order), false);
553 PHI_f = &*(N_face_bubble[0].data().begin());
554
555 int face_nodes[3] = {0, 1, 2};
557 face_nodes, face_order,
558 &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0), NULL, PHI_f_e, NULL,
559 nb_gauss_pts, 3, base_polynomials);
561 face_nodes, face_order,
562 &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0), NULL, PHI_f, NULL,
563 nb_gauss_pts, 3, base_polynomials);
564
565 // set shape functions into data structure
566 if (data.dataOnEntities[MBTRI].size() != 1) {
567 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
568 }
569 data.dataOnEntities[MBTRI][0].getN(base).resize(
570 nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_HDIV(face_order), false);
571 int col = 0;
572 for (int oo = 0; oo < face_order; oo++) {
573 for (int ee = 0; ee < 3; ee++) {
574 for (int dd = 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
575 dd < 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++, col++) {
576 for (int gg = 0; gg < nb_gauss_pts; gg++) {
577 data.dataOnEntities[MBTRI][0].getN(base)(gg, col) =
578 N_face_edge(0, ee)(gg, dd);
579 }
580 }
581 }
582 for (int dd = 3 * NBFACETRI_AINSWORTH_FACE_HDIV(oo);
583 dd < 3 * NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++, col++) {
584 for (int gg = 0; gg < nb_gauss_pts; gg++) {
585 data.dataOnEntities[MBTRI][0].getN(base)(gg, col) =
586 N_face_bubble[0](gg, dd);
587 }
588 }
589 }
590
592}
#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 594 of file TriPolynomialBase.cpp.

594 {
596
597 EntitiesFieldData &data = cTx->dAta;
598 // set shape functions into data structure
599 if (data.dataOnEntities[MBTRI].size() != 1) {
600 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
601 }
602 const FieldApproximationBase base = cTx->bAse;
603 if (base != DEMKOWICZ_JACOBI_BASE) {
604 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
605 "This should be used only with DEMKOWICZ_JACOBI_BASE "
606 "but base is %s",
608 }
609 int order = data.dataOnEntities[MBTRI][0].getOrder();
610 int nb_gauss_pts = pts.size2();
611 data.dataOnEntities[MBTRI][0].getN(base).resize(
612 nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
613 double *phi_f = &*data.dataOnEntities[MBTRI][0].getN(base).data().begin();
616 int face_nodes[3] = {0, 1, 2};
618 face_nodes, order,
619 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
620 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(), phi_f,
621 NULL, nb_gauss_pts, 3);
622
624}
#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 321 of file TriPolynomialBase.cpp.

321 {
323
324 switch (cTx->bAse) {
328 break;
331 break;
332 default:
333 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
334 }
335
337}
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)

◆ getValueL2AinsworthBase()

MoFEMErrorCode TriPolynomialBase::getValueL2AinsworthBase ( MatrixDouble pts)
private

Definition at line 339 of file TriPolynomialBase.cpp.

339 {
341
342 EntitiesFieldData &data = cTx->dAta;
343 const FieldApproximationBase base = cTx->bAse;
344 if (cTx->basePolynomialsType0 == NULL)
345 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
346 "Polynomial type not set");
347 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
348 double *diffL, const int dim) =
350
351 int nb_gauss_pts = pts.size2();
352
353 data.dataOnEntities[MBTRI][0].getN(base).resize(
354 nb_gauss_pts, NBFACETRI_L2(data.dataOnEntities[MBTRI][0].getOrder()),
355 false);
356 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(
357 nb_gauss_pts,
358 2 * NBFACETRI_L2(data.dataOnEntities[MBTRI][0].getOrder()), false);
359
361 data.dataOnEntities[MBTRI][0].getOrder(),
362 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
363 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
364 &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
365 &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
366 nb_gauss_pts, base_polynomials);
367
369}
#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 372 of file TriPolynomialBase.cpp.

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

10 {
12 *iface = const_cast<TriPolynomialBase *>(this);
14}
Calculate base functions on triangle.

Member Data Documentation

◆ cTx

EntPolynomialBaseCtx* MoFEM::TriPolynomialBase::cTx
private

Definition at line 28 of file TriPolynomialBase.hpp.

◆ diffN_face_bubble

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

Definition at line 41 of file TriPolynomialBase.hpp.

◆ diffN_face_edge

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

Definition at line 40 of file TriPolynomialBase.hpp.

◆ N_face_bubble

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

Definition at line 39 of file TriPolynomialBase.hpp.

◆ N_face_edge

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

Definition at line 38 of file TriPolynomialBase.hpp.


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