v0.13.1
EdgePolynomialBase.cpp
Go to the documentation of this file.
1/** \file EdgePolynomialBase.cpp
2\brief Implementation of Ainsworth-Cole H1 base on edge
3*/
4
5
6
7using namespace MoFEM;
8
10EdgePolynomialBase::query_interface(boost::typeindex::type_index type_index,
11 UnknownInterface **iface) const {
12 *iface = const_cast<EdgePolynomialBase *>(this);
13 return 0;
14}
15
18
21 boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
23
25
26 int nb_gauss_pts = pts.size2();
27 if (!nb_gauss_pts)
29
30 if (pts.size1() < 1)
31 SETERRQ(
32 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
33 "Wrong dimension of pts, should be at least 3 rows with coordinates");
34
35 const FieldApproximationBase base = cTx->bAse;
36 EntitiesFieldData &data = cTx->dAta;
37
39 if (cTx->copyNodeBase == LASTBASE) {
40 data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 2,
41 false);
43 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
44 &pts(0, 0), nb_gauss_pts);
45 } else
46 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
47 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
48
49 if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
50 (unsigned int)nb_gauss_pts)
51 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
52 "Base functions or nodes has wrong number of integration points "
53 "for base %s",
55
56 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(nb_gauss_pts, 2 * 1,
57 false);
58 for (auto gg = 0; gg != nb_gauss_pts; ++gg)
59 std::copy(Tools::diffShapeFunMBEDGE.begin(),
61 &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(gg, 0));
62 }
63
64 switch (cTx->sPace) {
65 case H1:
66 CHKERR getValueH1(pts);
67 break;
68 case HDIV:
70 break;
71 case HCURL:
73 break;
74 case L2:
75 CHKERR getValueL2(pts);
76 break;
77 default:
78 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
79 }
80
82}
83
86
87 switch (cTx->bAse) {
91 break;
94 break;
97 break;
98 default:
99 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
100 }
101
103}
104
108
109 EntitiesFieldData &data = cTx->dAta;
110 const std::string &field_name = cTx->fieldName;
111 const int nb_gauss_pts = pts.size2();
112
113 auto get_alpha = [field_name](auto &data) -> MatrixInt & {
114 auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
115 if (!ptr)
116 ptr.reset(new MatrixInt());
117 return *ptr;
118 };
119
120 auto get_base = [field_name](auto &data) -> MatrixDouble & {
121 auto &ptr = data.getBBNSharedPtr(field_name);
122 if (!ptr)
123 ptr.reset(new MatrixDouble());
124 return *ptr;
125 };
126
127 auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
128 auto &ptr = data.getBBDiffNSharedPtr(field_name);
129 if (!ptr)
130 ptr.reset(new MatrixDouble());
131 return *ptr;
132 };
133
134 auto &get_n = get_base(data.dataOnEntities[MBVERTEX][0]);
135 auto &get_diff_n = get_diff_base(data.dataOnEntities[MBVERTEX][0]);
136 get_n.resize(nb_gauss_pts, 2, false);
137 get_diff_n.resize(nb_gauss_pts, 2, false);
138 get_n.clear();
139 get_diff_n.clear();
140
141 if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
142 (unsigned int)nb_gauss_pts)
143 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
144 "Base functions or nodes has wrong number of integration points "
145 "for base %s",
147 auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
148
149 auto &vertex_alpha = get_alpha(data.dataOnEntities[MBVERTEX][0]);
150 vertex_alpha.resize(2, 2, false);
151 vertex_alpha(0, 0) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[0];
152 vertex_alpha(0, 1) = 0;
153 vertex_alpha(1, 0) = 0;
154 vertex_alpha(1, 1) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[1];
155
157 1, nb_gauss_pts, vertex_alpha.size1(), &vertex_alpha(0, 0), &lambda(0, 0),
158 Tools::diffShapeFunMBEDGE.data(), &get_n(0, 0), &get_diff_n(0, 0));
159 std::array<double, 2> f = {
160 boost::math::factorial<double>(
161 data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[0]),
162 boost::math::factorial<double>(
163 data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[1])};
164
165 for (int g = 0; g != nb_gauss_pts; ++g)
166 for (int n = 0; n != 2; ++n) {
167 get_n(g, n) *= f[n];
168 get_diff_n(g, n) *= f[n];
169 }
170
171 if (data.spacesOnEntities[MBEDGE].test(H1)) {
172 if (data.dataOnEntities[MBEDGE].size() != 1)
173 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
174 "Wrong size ent of ent data");
175
176 int order = data.dataOnEntities[MBEDGE][0].getOrder();
177 const int nb_dofs = NBEDGE_H1(order);
178
179 auto &get_n = get_base(data.dataOnEntities[MBEDGE][0]);
180 auto &get_diff_n = get_diff_base(data.dataOnEntities[MBEDGE][0]);
181 get_n.resize(nb_gauss_pts, nb_dofs, false);
182 get_diff_n.resize(nb_gauss_pts, nb_dofs, false);
183
184 if (nb_dofs) {
185 auto &edge_alpha = get_alpha(data.dataOnEntities[MBEDGE][0]);
186 edge_alpha.resize(nb_dofs, 2);
189 order, nb_gauss_pts, edge_alpha.size1(), &edge_alpha(0, 0),
190 &lambda(0, 0), Tools::diffShapeFunMBEDGE.data(), &get_n(0, 0),
191 &get_diff_n(0, 0));
192 }
193 } else {
194 get_n.resize(nb_gauss_pts, 0, false);
195 get_diff_n.resize(nb_gauss_pts, 0, false);
196 }
197
199}
200
203
204 EntitiesFieldData &data = cTx->dAta;
205 const FieldApproximationBase base = cTx->bAse;
206 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
207 double *diffL, const int dim) =
209
210 int nb_gauss_pts = pts.size2();
211
212 const int side_number = 0;
213 int sense = data.dataOnEntities[MBEDGE][side_number].getSense();
214 int order = data.dataOnEntities[MBEDGE][side_number].getOrder();
215 data.dataOnEntities[MBEDGE][side_number].getN(base).resize(
216 nb_gauss_pts, NBEDGE_H1(order), false);
217 data.dataOnEntities[MBEDGE][side_number].getDiffN(base).resize(
218 nb_gauss_pts, NBEDGE_H1(order), false);
219
220 data.dataOnEntities[MBEDGE][side_number].getN(base).clear();
221 data.dataOnEntities[MBEDGE][side_number].getDiffN(base).clear();
222
223 L.resize(NBEDGE_H1(order), false);
224 diffL.resize(NBEDGE_H1(order), false);
225
226 if (data.dataOnEntities[MBEDGE][side_number].getOrder() > 1) {
227
228 double diff_s = 2.; // s = s(xi), ds/dxi = 2., because change of basis
229 for (int gg = 0; gg != nb_gauss_pts; gg++) {
230
231 double s = 2 * pts(0, gg) - 1; // makes form -1..1
232 if (!sense) {
233 s *= -1;
234 diff_s *= -1;
235 }
236
237 // calculate Legendre polynomials at integration points
238 CHKERR base_polynomials(NBEDGE_H1(order) - 1, s, &diff_s,
239 &*L.data().begin(), &*diffL.data().begin(), 1);
240
241 for (unsigned int pp = 0;
242 pp < data.dataOnEntities[MBEDGE][side_number].getN(base).size2();
243 pp++) {
244
245 // Calculate edge shape functions N0*N1*L(p), where N0 and N1 are nodal
246 // shape functions
247 double v = data.dataOnEntities[MBVERTEX][0].getN(base)(gg, 0) *
248 data.dataOnEntities[MBVERTEX][0].getN(base)(gg, 1);
249 data.dataOnEntities[MBEDGE][side_number].getN(base)(gg, pp) = v * L(pp);
250
251 // Calculate derivative edge shape functions
252 // dN/dksi = dN0/dxi*N1*L + N0*dN1/ksi*L + N0*N1*dL/dxi
253 data.dataOnEntities[MBEDGE][side_number].getDiffN(base)(gg, pp) =
254 ((+1.) * data.dataOnEntities[MBVERTEX][0].getN(base)(gg, 1) +
255 data.dataOnEntities[MBVERTEX][0].getN(base)(gg, 0) * (-1.)) *
256 L(pp) +
257 v * diffL(pp);
258 }
259 }
260 }
261
263}
264
267
268 EntitiesFieldData &data = cTx->dAta;
269 const FieldApproximationBase base = cTx->bAse;
270
271 int nb_gauss_pts = pts.size2();
272
273 const int side_number = 0;
274 int order = data.dataOnEntities[MBEDGE][side_number].getOrder();
275
276 data.dataOnEntities[MBEDGE][side_number].getN(base).resize(
277 nb_gauss_pts, NBEDGE_H1(order), false);
278 data.dataOnEntities[MBEDGE][side_number].getDiffN(base).resize(
279 nb_gauss_pts, NBEDGE_H1(order), false);
280
281 data.dataOnEntities[MBEDGE][side_number].getN(base).clear();
282 data.dataOnEntities[MBEDGE][side_number].getDiffN(base).clear();
283
284 double diff_shape[] = {-1, 1};
285 MatrixDouble shape(nb_gauss_pts, 2);
286
287 if (data.dataOnEntities[MBEDGE][side_number].getOrder() > 1) {
288 for (int gg = 0; gg != nb_gauss_pts; gg++) {
289 const double mu0 = 1.0 - pts(0, gg); // pts ranges over [0, 1]
290 const double mu1 = pts(0, gg);
291 shape(gg, 0) = mu0;
292 shape(gg, 1) = mu1;
293 }
294
296 order, &*shape.data().begin(), diff_shape,
297 &*data.dataOnEntities[MBEDGE][side_number].getN(base).data().begin(),
298 &*data.dataOnEntities[MBEDGE][side_number]
299 .getDiffN(base)
300 .data()
301 .begin(),
302 nb_gauss_pts);
303 }
304
306}
307
310
311 switch (cTx->bAse) {
314 break;
315 default:
316 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
317 }
318
320}
321
324
325 EntitiesFieldData &data = cTx->dAta;
326 const FieldApproximationBase base = cTx->bAse;
327
328 int nb_gauss_pts = pts.size2();
329
330 const int side_number = 0;
331 int order = data.dataOnEntities[MBEDGE][side_number].getOrder();
332
333 data.dataOnEntities[MBEDGE][side_number].getN(base).resize(
334 nb_gauss_pts, NBEDGE_L2(order), false);
335 data.dataOnEntities[MBEDGE][side_number].getDiffN(base).resize(
336 nb_gauss_pts, NBEDGE_L2(order), false);
337
338 data.dataOnEntities[MBEDGE][side_number].getN(base).clear();
339
340 double diff_n[] = {-1, 1};
341 MatrixDouble shape(nb_gauss_pts, 2);
342 if (data.dataOnEntities[MBEDGE][side_number].getOrder() > 1) {
343 for (int gg = 0; gg != nb_gauss_pts; gg++) {
344 const double mu0 = 1.0 - pts(0, gg); // pts ranges over [0, 1]
345 const double mu1 = pts(0, gg);
346 shape(gg, 0) = mu0;
347 shape(gg, 1) = mu1;
348 }
349
351 order, &*shape.data().begin(), diff_n,
352 &*data.dataOnEntities[MBEDGE][side_number].getN(base).data().begin(),
353 &*data.dataOnEntities[MBEDGE][side_number]
354 .getDiffN(base)
355 .data()
356 .begin(),
357 nb_gauss_pts);
358 }
359
361}
362
364 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
365 "Make no sense, unless problem is 2d (2d not implemented yet)");
366}
367
371
372 EntitiesFieldData &data = cTx->dAta;
373 const FieldApproximationBase base = cTx->bAse;
374
375 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
376 double *diffL, const int dim) =
378
379 int nb_gauss_pts = pts.size2();
380 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
381
382 if (data.dataOnEntities[MBEDGE].size() != 1)
383 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
384
385 int sense = data.dataOnEntities[MBEDGE][0].getSense();
386 int order = data.dataOnEntities[MBEDGE][0].getOrder();
387 int nb_dofs =
388 NBEDGE_AINSWORTH_HCURL(data.dataOnEntities[MBEDGE][0].getOrder());
389 data.dataOnEntities[MBEDGE][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
390 false);
391 data.dataOnEntities[MBEDGE][0].getDiffN(base).resize(nb_gauss_pts, 0,
392 false);
393
395 sense, order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
396 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
397 &*data.dataOnEntities[MBEDGE][0].getN(base).data().begin(), nullptr,
398 nb_gauss_pts, base_polynomials);
399
400 } else {
401 data.dataOnEntities[MBEDGE][0].getN(base).resize(nb_gauss_pts, 0, false);
402 data.dataOnEntities[MBEDGE][0].getDiffN(base).resize(nb_gauss_pts, 0,
403 false);
404 }
405
407}
408
412
413 EntitiesFieldData &data = cTx->dAta;
414 const FieldApproximationBase base = cTx->bAse;
415
416 int nb_gauss_pts = pts.size2();
417 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
418
419 if (data.dataOnEntities[MBEDGE].size() != 1)
420 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
421 "No data structure to store base functions");
422
423 int sense = data.dataOnEntities[MBEDGE][0].getSense();
424 int order = data.dataOnEntities[MBEDGE][0].getOrder();
425 int nb_dofs =
426 NBEDGE_DEMKOWICZ_HCURL(data.dataOnEntities[MBEDGE][0].getOrder());
427 data.dataOnEntities[MBEDGE][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
428 false);
429 data.dataOnEntities[MBEDGE][0].getDiffN(base).resize(nb_gauss_pts, 0,
430 false);
432 sense, order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
433 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
434 &*data.dataOnEntities[MBEDGE][0].getN(base).data().begin(), NULL,
435 nb_gauss_pts);
436 } else {
437
438 data.dataOnEntities[MBEDGE][0].getN(base).resize(nb_gauss_pts, 0, false);
439 data.dataOnEntities[MBEDGE][0].getDiffN(base).resize(nb_gauss_pts, 0,
440 false);
441 }
442
444}
445
448
449 switch (cTx->bAse) {
453 break;
456 break;
457 default:
458 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
459 }
460
462}
static Index< 'p', 3 > p
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
@ NOBASE
Definition: definitions.h:59
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ 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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
const int dim
PetscErrorCode ShapeMBEDGE(double *N, const double *G_X, int DIM)
Definition: fem_tools.c:761
FTensor::Index< 'n', SPACE_DIM > n
constexpr double lambda
#define NBEDGE_DEMKOWICZ_HCURL(P)
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
#define NBEDGE_L2(P)
Number of base functions on edge fro L2 space.
#define NBEDGE_AINSWORTH_HCURL(P)
const double v
phase velocity of light in medium (cm/ns)
auto f
Definition: HenckyOps.hpp:5
MoFEMErrorCode L2_ShapeFunctions_ONSEGMENT(int p, double *N, double *diffN, double *funN, double *funDiffN, int nb_integration_pts)
MoFEMErrorCode H1_BubbleShapeFunctions_ONSEGMENT(int p, double *N, double *diffN, double *bubbleN, double *diff_bubbleN, int nb_integration_pts)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasMatrix< int > MatrixInt
Definition: Types.hpp:76
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
MoFEMErrorCode Hcurl_Demkowicz_EdgeBaseFunctions_MBEDGE(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 edge.
Definition: Hcurl.cpp:2144
MoFEMErrorCode Hcurl_Ainsworth_EdgeBaseFunctions_MBTET_ON_EDGE(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 edge.
Definition: Hcurl.cpp:175
constexpr auto field_name
constexpr double g
static MoFEMErrorCode baseFunctionsEdge(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 generateIndicesEdgeEdge(const int N, int *alpha)
Calculate base functions on tetrahedral.
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
MoFEMErrorCode getValueH1DemkowiczBase(MatrixDouble &pts)
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)
MoFEMErrorCode getValueL2(MatrixDouble &pts)
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
EntPolynomialBaseCtx * cTx
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)
MoFEMErrorCode getValueH1(MatrixDouble &pts)
MoFEMErrorCode getValueL2DemkowiczBase(MatrixDouble &pts)
Class used to pass element data to calculate base functions on tet,triangle,edge.
const FieldApproximationBase copyNodeBase
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
const FieldApproximationBase bAse
data structure for finite element entity
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
static constexpr std::array< double, 2 > diffShapeFunMBEDGE
Definition: Tools.hpp:68
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.