v0.14.0
Loading...
Searching...
No Matches
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 boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
20
22
23 int nb_gauss_pts = pts.size2();
24 if (!nb_gauss_pts)
26
27 if (pts.size1() < 1)
28 SETERRQ(
29 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
30 "Wrong dimension of pts, should be at least 3 rows with coordinates");
31
32 const FieldApproximationBase base = cTx->bAse;
33 EntitiesFieldData &data = cTx->dAta;
34
36 if (cTx->copyNodeBase == LASTBASE) {
37 data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 2,
38 false);
40 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
41 &pts(0, 0), nb_gauss_pts);
42 } else
43 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
44 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
45
46 if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
47 (unsigned int)nb_gauss_pts)
48 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
49 "Base functions or nodes has wrong number of integration points "
50 "for base %s",
52
53 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(nb_gauss_pts, 2 * 1,
54 false);
55 for (auto gg = 0; gg != nb_gauss_pts; ++gg)
56 std::copy(Tools::diffShapeFunMBEDGE.begin(),
58 &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(gg, 0));
59 }
60
61 switch (cTx->sPace) {
62 case H1:
63 CHKERR getValueH1(pts);
64 break;
65 case HDIV:
67 break;
68 case HCURL:
70 break;
71 case L2:
72 CHKERR getValueL2(pts);
73 break;
74 default:
75 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
76 }
77
79}
80
83
84 switch (cTx->bAse) {
88 break;
91 break;
94 break;
95 default:
96 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
97 }
98
100}
101
105
106 EntitiesFieldData &data = cTx->dAta;
107 const std::string &field_name = cTx->fieldName;
108 const int nb_gauss_pts = pts.size2();
109
110 auto get_alpha = [field_name](auto &data) -> MatrixInt & {
111 auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
112 if (!ptr)
113 ptr.reset(new MatrixInt());
114 return *ptr;
115 };
116
117 auto get_base = [field_name](auto &data) -> MatrixDouble & {
118 auto &ptr = data.getBBNSharedPtr(field_name);
119 if (!ptr)
120 ptr.reset(new MatrixDouble());
121 return *ptr;
122 };
123
124 auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
125 auto &ptr = data.getBBDiffNSharedPtr(field_name);
126 if (!ptr)
127 ptr.reset(new MatrixDouble());
128 return *ptr;
129 };
130
131 auto &get_n = get_base(data.dataOnEntities[MBVERTEX][0]);
132 auto &get_diff_n = get_diff_base(data.dataOnEntities[MBVERTEX][0]);
133 get_n.resize(nb_gauss_pts, 2, false);
134 get_diff_n.resize(nb_gauss_pts, 2, false);
135 get_n.clear();
136 get_diff_n.clear();
137
138 if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
139 (unsigned int)nb_gauss_pts)
140 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
141 "Base functions or nodes has wrong number of integration points "
142 "for base %s",
144 auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
145
146 auto &vertex_alpha = get_alpha(data.dataOnEntities[MBVERTEX][0]);
147 vertex_alpha.resize(2, 2, false);
148 vertex_alpha(0, 0) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[0];
149 vertex_alpha(0, 1) = 0;
150 vertex_alpha(1, 0) = 0;
151 vertex_alpha(1, 1) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[1];
152
154 1, nb_gauss_pts, vertex_alpha.size1(), &vertex_alpha(0, 0), &lambda(0, 0),
155 Tools::diffShapeFunMBEDGE.data(), &get_n(0, 0), &get_diff_n(0, 0));
156 std::array<double, 2> f = {
157 boost::math::factorial<double>(
158 data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[0]),
159 boost::math::factorial<double>(
160 data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[1])};
161
162 for (int g = 0; g != nb_gauss_pts; ++g)
163 for (int n = 0; n != 2; ++n) {
164 get_n(g, n) *= f[n];
165 get_diff_n(g, n) *= f[n];
166 }
167
168 if (data.spacesOnEntities[MBEDGE].test(H1)) {
169 if (data.dataOnEntities[MBEDGE].size() != 1)
170 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
171 "Wrong size ent of ent data");
172
173 int order = data.dataOnEntities[MBEDGE][0].getOrder();
174 const int nb_dofs = NBEDGE_H1(order);
175
176 auto &get_n = get_base(data.dataOnEntities[MBEDGE][0]);
177 auto &get_diff_n = get_diff_base(data.dataOnEntities[MBEDGE][0]);
178 get_n.resize(nb_gauss_pts, nb_dofs, false);
179 get_diff_n.resize(nb_gauss_pts, nb_dofs, false);
180
181 if (nb_dofs) {
182 auto &edge_alpha = get_alpha(data.dataOnEntities[MBEDGE][0]);
183 edge_alpha.resize(nb_dofs, 2);
186 order, nb_gauss_pts, edge_alpha.size1(), &edge_alpha(0, 0),
187 &lambda(0, 0), Tools::diffShapeFunMBEDGE.data(), &get_n(0, 0),
188 &get_diff_n(0, 0));
189 }
190 } else {
191 get_n.resize(nb_gauss_pts, 0, false);
192 get_diff_n.resize(nb_gauss_pts, 0, false);
193 }
194
196}
197
200
201 EntitiesFieldData &data = cTx->dAta;
202 const FieldApproximationBase base = cTx->bAse;
203 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
204 double *diffL, const int dim) =
206
207 int nb_gauss_pts = pts.size2();
208
209 const int side_number = 0;
210 int sense = data.dataOnEntities[MBEDGE][side_number].getSense();
211 int order = data.dataOnEntities[MBEDGE][side_number].getOrder();
212 data.dataOnEntities[MBEDGE][side_number].getN(base).resize(
213 nb_gauss_pts, NBEDGE_H1(order), false);
214 data.dataOnEntities[MBEDGE][side_number].getDiffN(base).resize(
215 nb_gauss_pts, NBEDGE_H1(order), false);
216
217 data.dataOnEntities[MBEDGE][side_number].getN(base).clear();
218 data.dataOnEntities[MBEDGE][side_number].getDiffN(base).clear();
219
220 L.resize(NBEDGE_H1(order), false);
221 diffL.resize(NBEDGE_H1(order), false);
222
223 if (data.dataOnEntities[MBEDGE][side_number].getOrder() > 1) {
224
225 double diff_s = 2.; // s = s(xi), ds/dxi = 2., because change of basis
226 for (int gg = 0; gg != nb_gauss_pts; gg++) {
227
228 double s = 2 * pts(0, gg) - 1; // makes form -1..1
229 if (!sense) {
230 s *= -1;
231 diff_s *= -1;
232 }
233
234 // calculate Legendre polynomials at integration points
235 CHKERR base_polynomials(NBEDGE_H1(order) - 1, s, &diff_s,
236 &*L.data().begin(), &*diffL.data().begin(), 1);
237
238 for (unsigned int pp = 0;
239 pp < data.dataOnEntities[MBEDGE][side_number].getN(base).size2();
240 pp++) {
241
242 // Calculate edge shape functions N0*N1*L(p), where N0 and N1 are nodal
243 // shape functions
244 double v = data.dataOnEntities[MBVERTEX][0].getN(base)(gg, 0) *
245 data.dataOnEntities[MBVERTEX][0].getN(base)(gg, 1);
246 data.dataOnEntities[MBEDGE][side_number].getN(base)(gg, pp) = v * L(pp);
247
248 // Calculate derivative edge shape functions
249 // dN/dksi = dN0/dxi*N1*L + N0*dN1/ksi*L + N0*N1*dL/dxi
250 data.dataOnEntities[MBEDGE][side_number].getDiffN(base)(gg, pp) =
251 ((+1.) * data.dataOnEntities[MBVERTEX][0].getN(base)(gg, 1) +
252 data.dataOnEntities[MBVERTEX][0].getN(base)(gg, 0) * (-1.)) *
253 L(pp) +
254 v * diffL(pp);
255 }
256 }
257 }
258
260}
261
264
265 EntitiesFieldData &data = cTx->dAta;
266 const FieldApproximationBase base = cTx->bAse;
267
268 int nb_gauss_pts = pts.size2();
269
270 const int side_number = 0;
271 int order = data.dataOnEntities[MBEDGE][side_number].getOrder();
272
273 data.dataOnEntities[MBEDGE][side_number].getN(base).resize(
274 nb_gauss_pts, NBEDGE_H1(order), false);
275 data.dataOnEntities[MBEDGE][side_number].getDiffN(base).resize(
276 nb_gauss_pts, NBEDGE_H1(order), false);
277
278 data.dataOnEntities[MBEDGE][side_number].getN(base).clear();
279 data.dataOnEntities[MBEDGE][side_number].getDiffN(base).clear();
280
281 double diff_shape[] = {-1, 1};
282 MatrixDouble shape(nb_gauss_pts, 2);
283
284 if (data.dataOnEntities[MBEDGE][side_number].getOrder() > 1) {
285 for (int gg = 0; gg != nb_gauss_pts; gg++) {
286 const double mu0 = 1.0 - pts(0, gg); // pts ranges over [0, 1]
287 const double mu1 = pts(0, gg);
288 shape(gg, 0) = mu0;
289 shape(gg, 1) = mu1;
290 }
291
293 order, &*shape.data().begin(), diff_shape,
294 &*data.dataOnEntities[MBEDGE][side_number].getN(base).data().begin(),
295 &*data.dataOnEntities[MBEDGE][side_number]
296 .getDiffN(base)
297 .data()
298 .begin(),
299 nb_gauss_pts);
300 }
301
303}
304
307
308 switch (cTx->bAse) {
311 break;
312 default:
313 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
314 }
315
317}
318
321
322 EntitiesFieldData &data = cTx->dAta;
323 const FieldApproximationBase base = cTx->bAse;
324
325 int nb_gauss_pts = pts.size2();
326
327 const int side_number = 0;
328 int order = data.dataOnEntities[MBEDGE][side_number].getOrder();
329
330 data.dataOnEntities[MBEDGE][side_number].getN(base).resize(
331 nb_gauss_pts, NBEDGE_L2(order), false);
332 data.dataOnEntities[MBEDGE][side_number].getDiffN(base).resize(
333 nb_gauss_pts, NBEDGE_L2(order), false);
334
335 data.dataOnEntities[MBEDGE][side_number].getN(base).clear();
336
337 double diff_n[] = {-1, 1};
338 MatrixDouble shape(nb_gauss_pts, 2);
339 if (data.dataOnEntities[MBEDGE][side_number].getOrder() > 1) {
340 for (int gg = 0; gg != nb_gauss_pts; gg++) {
341 const double mu0 = 1.0 - pts(0, gg); // pts ranges over [0, 1]
342 const double mu1 = pts(0, gg);
343 shape(gg, 0) = mu0;
344 shape(gg, 1) = mu1;
345 }
346
348 order, &*shape.data().begin(), diff_n,
349 &*data.dataOnEntities[MBEDGE][side_number].getN(base).data().begin(),
350 &*data.dataOnEntities[MBEDGE][side_number]
351 .getDiffN(base)
352 .data()
353 .begin(),
354 nb_gauss_pts);
355 }
356
358}
359
361 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
362 "Make no sense, unless problem is 2d (2d not implemented yet)");
363}
364
368
369 EntitiesFieldData &data = cTx->dAta;
370 const FieldApproximationBase base = cTx->bAse;
371
372 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
373 double *diffL, const int dim) =
375
376 int nb_gauss_pts = pts.size2();
377 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
378
379 if (data.dataOnEntities[MBEDGE].size() != 1)
380 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
381
382 int sense = data.dataOnEntities[MBEDGE][0].getSense();
383 int order = data.dataOnEntities[MBEDGE][0].getOrder();
384 int nb_dofs =
385 NBEDGE_AINSWORTH_HCURL(data.dataOnEntities[MBEDGE][0].getOrder());
386 data.dataOnEntities[MBEDGE][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
387 false);
388 data.dataOnEntities[MBEDGE][0].getDiffN(base).resize(nb_gauss_pts, 0,
389 false);
390
392 sense, order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
393 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
394 &*data.dataOnEntities[MBEDGE][0].getN(base).data().begin(), nullptr,
395 nb_gauss_pts, base_polynomials);
396
397 } else {
398 data.dataOnEntities[MBEDGE][0].getN(base).resize(nb_gauss_pts, 0, false);
399 data.dataOnEntities[MBEDGE][0].getDiffN(base).resize(nb_gauss_pts, 0,
400 false);
401 }
402
404}
405
409
410 EntitiesFieldData &data = cTx->dAta;
411 const FieldApproximationBase base = cTx->bAse;
412
413 int nb_gauss_pts = pts.size2();
414 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
415
416 if (data.dataOnEntities[MBEDGE].size() != 1)
417 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
418 "No data structure to store base functions");
419
420 int sense = data.dataOnEntities[MBEDGE][0].getSense();
421 int order = data.dataOnEntities[MBEDGE][0].getOrder();
422 int nb_dofs =
423 NBEDGE_DEMKOWICZ_HCURL(data.dataOnEntities[MBEDGE][0].getOrder());
424 data.dataOnEntities[MBEDGE][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
425 false);
426 data.dataOnEntities[MBEDGE][0].getDiffN(base).resize(nb_gauss_pts, 0,
427 false);
429 sense, order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
430 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
431 &*data.dataOnEntities[MBEDGE][0].getN(base).data().begin(), NULL,
432 nb_gauss_pts);
433 } else {
434
435 data.dataOnEntities[MBEDGE][0].getN(base).resize(nb_gauss_pts, 0, false);
436 data.dataOnEntities[MBEDGE][0].getDiffN(base).resize(nb_gauss_pts, 0,
437 false);
438 }
439
441}
442
445
446 switch (cTx->bAse) {
450 break;
453 break;
454 default:
455 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
456 }
457
459}
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
constexpr int order
const int dim
PetscErrorCode ShapeMBEDGE(double *N, const double *G_X, int DIM)
Definition: fem_tools.c:761
FTensor::Index< 'n', SPACE_DIM > n
#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)
static double lambda
const double v
phase velocity of light in medium (cm/ns)
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: Common.hpp:10
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 query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
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.