v0.14.0
Loading...
Searching...
No Matches
DataOperators.cpp
Go to the documentation of this file.
1/** file DataOperators.cpp
2
3 \brief implementation of Data Operators for Forces and Sources
4
5*/
6
7
8
9#ifdef __cplusplus
10extern "C" {
11#endif
12#include <cblas.h>
13#include <lapack_wrap.h>
14#include <gm_rule.h>
15#ifdef __cplusplus
16}
17#endif
18
19namespace MoFEM {
20
22 :
23
24 sYmm(symm),
25
26 doEntities{true, true, true, true, true, true,
27 true, true, true, true, true, true},
28
29 doVertices(doEntities[MBVERTEX]), doEdges(doEntities[MBEDGE]),
30 doQuads(doEntities[MBQUAD]), doTris(doEntities[MBTRI]),
31 doTets(doEntities[MBTET]), doPrisms(doEntities[MBPRISM]) {
32
33 /// This not yet implemented, switch off.
34 doEntities[MBPOLYGON] = false;
35 doEntities[MBPYRAMID] = false;
36 doEntities[MBKNIFE] = false;
37 doEntities[MBPOLYHEDRON] = false;
38}
39
40template <bool Symm>
42 EntitiesFieldData &col_data) {
44
45 auto do_col_entity =
46 [&](boost::ptr_vector<EntitiesFieldData::EntData> &row_ent_data,
47 const int ss, const EntityType row_type, const EntityType low_type,
48 const EntityType hi_type) {
50 for (EntityType col_type = low_type; col_type != hi_type; ++col_type) {
51 auto &col_ent_data = col_data.dataOnEntities[col_type];
52 for (size_t SS = 0; SS != col_ent_data.size(); SS++) {
53 if (col_ent_data[SS].getFieldData().size())
54 CHKERR doWork(ss, SS, row_type, col_type, row_ent_data[ss],
55 col_ent_data[SS]);
56 }
57 }
59 };
60
61 auto do_row_entity = [&](const EntityType type) {
63 auto &row_ent_data = row_data.dataOnEntities[type];
64 for (size_t ss = 0; ss != row_ent_data.size(); ++ss) {
65 if constexpr (!Symm)
66 CHKERR do_col_entity(row_ent_data, ss, type, MBVERTEX, type);
67 size_t SS = 0;
68 if constexpr (Symm)
69 SS = ss;
70 for (; SS < col_data.dataOnEntities[type].size(); ++SS) {
71 CHKERR doWork(ss, SS, type, type, row_ent_data[ss],
72 col_data.dataOnEntities[type][SS]);
73 }
74 CHKERR do_col_entity(row_ent_data, ss, type,
75 static_cast<EntityType>(type + 1), MBMAXTYPE);
76 }
78 };
79
80 for (EntityType row_type = MBVERTEX; row_type != MBMAXTYPE; ++row_type) {
81 if (doEntities[row_type]) {
82 CHKERR do_row_entity(row_type);
83 }
84 }
85
86
88}
89
91 EntitiesFieldData &col_data) {
92 if (getSymm())
93 return opLhs<true>(row_data, col_data);
94 else
95 return opLhs<false>(row_data, col_data);
96}
97
98template <bool ErrorIfNoBase>
101 const std::array<bool, MBMAXTYPE> &do_entities) {
103
104 auto do_entity = [&](auto type) {
106
107 auto &ent_data = data.dataOnEntities[type];
108 const size_t size = ent_data.size();
109 for (size_t ss = 0; ss != size; ++ss) {
110
111 auto &side_data = ent_data[ss];
112
113 if (ErrorIfNoBase) {
114 if (side_data.getFieldData().size() &&
115 (side_data.getBase() == NOBASE ||
116 side_data.getBase() == LASTBASE)) {
117 for (VectorDofs::iterator it = side_data.getFieldDofs().begin();
118 it != side_data.getFieldDofs().end(); it++)
119 if ((*it) && (*it)->getActive())
120 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "No base on");
121 }
122 }
123
124 CHKERR doWork(ss, type, side_data);
125 }
126
128 };
129
130 for (EntityType row_type = MBVERTEX; row_type != MBMAXTYPE; ++row_type) {
131 if (do_entities[row_type]) {
132 CHKERR do_entity(row_type);
133 }
134 }
135
136
138}
139
141 const bool error_if_no_base) {
142 if (error_if_no_base)
143 return opRhs<true>(data, doEntities);
144 else
145 return opRhs<false>(data, doEntities);
146}
147
148template <>
150 MatrixDouble &jac_data, VectorDouble &det_data,
151 MatrixDouble &inv_jac_data) {
153 auto A = getFTensor2FromMat<3, 3>(jac_data);
154 int nb_gauss_pts = jac_data.size2();
155 det_data.resize(nb_gauss_pts, false);
156 inv_jac_data.resize(3, nb_gauss_pts, false);
157 auto det = getFTensor0FromVec(det_data);
158 auto I = getFTensor2FromMat<3, 3>(inv_jac_data);
159 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
161 CHKERR invertTensor3by3(A, det, I);
162 ++A;
163 ++det;
164 ++I;
165 }
167}
168
172
173 auto transform_base = [&](MatrixDouble &diff_n) {
175
176 if (diff_n.data().size()) {
177 const int nb_base_functions = diff_n.size2() / 3;
178 const int nb_gauss_pts = diff_n.size1();
179 diffNinvJac.resize(diff_n.size1(), diff_n.size2(), false);
180
181 double *t_diff_n_ptr = &*diff_n.data().begin();
183 t_diff_n_ptr, &t_diff_n_ptr[1], &t_diff_n_ptr[2]);
184 double *t_inv_n_ptr = &*diffNinvJac.data().begin();
186 t_inv_n_ptr, &t_inv_n_ptr[1], &t_inv_n_ptr[2]);
187
188 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
189 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
190 t_inv_diff_n(i) = t_diff_n(j) * tInvJac(j, i);
191 ++t_diff_n;
192 ++t_inv_diff_n;
193 }
194 }
195 diff_n.swap(diffNinvJac);
196 }
197
199 };
200
201 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
202 const auto base = static_cast<FieldApproximationBase>(b);
203 CHKERR transform_base(data.getDiffN(base));
204 }
205
206 switch (type) {
207 case MBVERTEX:
208 for (auto &m : data.getBBDiffNMap())
209 if (m.second)
210 CHKERR transform_base(*(m.second));
211 break;
212 default:
213 for (auto &ptr : data.getBBDiffNByOrderArray())
214 if (ptr)
215 CHKERR transform_base(*ptr);
216 }
217
219}
220
225
226 if (type == MBVERTEX)
228
229 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
230
231 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
232
233 const unsigned int nb_gauss_pts = data.getDiffN(base).size1();
234 const unsigned int nb_base_functions = data.getDiffN(base).size2() / 9;
235 if (!nb_base_functions)
236 continue;
237
238 diffHdivInvJac.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
239
240 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
241 double *t_inv_diff_n_ptr = &*diffHdivInvJac.data().begin();
243 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[HVEC0_1],
244 &t_inv_diff_n_ptr[HVEC0_2],
245
246 &t_inv_diff_n_ptr[HVEC1_0], &t_inv_diff_n_ptr[HVEC1_1],
247 &t_inv_diff_n_ptr[HVEC1_2],
248
249 &t_inv_diff_n_ptr[HVEC2_0], &t_inv_diff_n_ptr[HVEC2_1],
250 &t_inv_diff_n_ptr[HVEC2_2]);
251
252 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
253 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
254 t_inv_diff_n(k, i) = t_diff_n(k, j) * tInvJac(j, i);
255 ++t_diff_n;
256 ++t_inv_diff_n;
257 }
258 }
259
260 data.getDiffN(base).swap(diffHdivInvJac);
261 }
262
264}
265
267 int side, EntityType type, EntitiesFieldData::EntData &data) {
269
270 if (CN::Dimension(type) > 1) {
271
272 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
273
274 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
275
276 const unsigned int nb_base_functions = data.getN(base).size2() / 3;
277 if (!nb_base_functions)
278 continue;
279
280 const unsigned int nb_gauss_pts = data.getN(base).size1();
281 double const a = 1. / vOlume;
282
283 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
284 if (data.getN(base).size2() > 0) {
285 auto t_n = data.getFTensor1N<3>(base);
286 double *t_transformed_n_ptr = &*piolaN.data().begin();
288 t_transformed_n_ptr, // HVEC0
289 &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
290 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
291 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
292 t_transformed_n(i) = a * (tJac(i, k) * t_n(k));
293 ++t_n;
294 ++t_transformed_n;
295 }
296 }
297 data.getN(base).swap(piolaN);
298 }
299
300 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
301 if (data.getDiffN(base).size2() > 0) {
302 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
303 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
305 t_transformed_diff_n(t_transformed_diff_n_ptr,
306 &t_transformed_diff_n_ptr[HVEC0_1],
307 &t_transformed_diff_n_ptr[HVEC0_2],
308 &t_transformed_diff_n_ptr[HVEC1_0],
309 &t_transformed_diff_n_ptr[HVEC1_1],
310 &t_transformed_diff_n_ptr[HVEC1_2],
311 &t_transformed_diff_n_ptr[HVEC2_0],
312 &t_transformed_diff_n_ptr[HVEC2_1],
313 &t_transformed_diff_n_ptr[HVEC2_2]);
314 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
315 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
316 t_transformed_diff_n(i, k) = a * tJac(i, j) * t_diff_n(j, k);
317 ++t_diff_n;
318 ++t_transformed_diff_n;
319 }
320 }
321 data.getDiffN(base).swap(piolaDiffN);
322 }
323 }
324 }
325
327}
328
333
334 if (type == MBVERTEX)
336
337 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
338
339 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
340
341 const unsigned int nb_base_functions = data.getN(base).size2() / 3;
342 if (!nb_base_functions)
343 continue;
344
345 const unsigned int nb_gauss_pts = data.getN(base).size1();
346 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
347 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
348
349 auto t_n = data.getFTensor1N<3>(base);
350 double *t_transformed_n_ptr = &*piolaN.data().begin();
352 t_transformed_n_ptr, &t_transformed_n_ptr[HVEC1],
353 &t_transformed_n_ptr[HVEC2]);
354 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
355 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
356 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_transformed_diff_n(
357 t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[HVEC0_1],
358 &t_transformed_diff_n_ptr[HVEC0_2], &t_transformed_diff_n_ptr[HVEC1_0],
359 &t_transformed_diff_n_ptr[HVEC1_1], &t_transformed_diff_n_ptr[HVEC1_2],
360 &t_transformed_diff_n_ptr[HVEC2_0], &t_transformed_diff_n_ptr[HVEC2_1],
361 &t_transformed_diff_n_ptr[HVEC2_2]);
362
363 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
364 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
365 t_transformed_n(i) = tInvJac(k, i) * t_n(k);
366 t_transformed_diff_n(i, k) = tInvJac(j, i) * t_diff_n(j, k);
367 ++t_n;
368 ++t_transformed_n;
369 ++t_diff_n;
370 ++t_transformed_diff_n;
371 }
372 }
373 data.getN(base).swap(piolaN);
374 data.getDiffN(base).swap(piolaDiffN);
375 }
376
377 // data.getBase() = base;
378
380}
381
386
387 if (data.getFieldData().size() == 0)
389 const int valid_edges3[] = {1, 1, 1, 0, 0, 0, 0, 0, 0};
390 const int valid_faces3[] = {0, 0, 0, 1, 0, 0, 0, 0, 0};
391 const int valid_edges4[] = {0, 0, 0, 0, 0, 0, 1, 1, 1};
392 const int valid_faces4[] = {0, 0, 0, 0, 1, 0, 0, 0, 0};
393
394 if (type == MBEDGE) {
395 if (!valid_edges3[side] || valid_edges4[side])
397 } else if (type == MBTRI) {
398 if (!valid_faces3[side] || valid_faces4[side])
400 }
401
402 switch (type) {
403 case MBVERTEX: {
404 for (unsigned int gg = 0; gg < data.getN().size1(); ++gg) {
405 for (int dd = 0; dd < 3; dd++) {
406 cOords_at_GaussPtF3(gg, dd) =
407 cblas_ddot(3, &data.getN(gg)[0], 1, &data.getFieldData()[dd], 3);
408 tAngent1_at_GaussPtF3(gg, dd) = cblas_ddot(
409 3, &data.getDiffN()(gg, 0), 2, &data.getFieldData()[dd], 3);
410 tAngent2_at_GaussPtF3(gg, dd) = cblas_ddot(
411 3, &data.getDiffN()(gg, 1), 2, &data.getFieldData()[dd], 3);
412 cOords_at_GaussPtF4(gg, dd) = cblas_ddot(
413 3, &data.getN(gg)[0], 1, &data.getFieldData()[9 + dd], 3);
414 tAngent1_at_GaussPtF4(gg, dd) = cblas_ddot(
415 3, &data.getDiffN()(gg, 6 + 0), 2, &data.getFieldData()[9 + dd], 3);
416 tAngent2_at_GaussPtF4(gg, dd) = cblas_ddot(
417 3, &data.getDiffN()(gg, 6 + 1), 2, &data.getFieldData()[9 + dd], 3);
418 }
419 }
420 } break;
421 case MBEDGE:
422 case MBTRI: {
423 if (2 * data.getN().size2() != data.getDiffN().size2()) {
424 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
425 }
426 unsigned int nb_dofs = data.getFieldData().size();
427 if (nb_dofs % 3 != 0) {
428 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
429 }
430 if (nb_dofs > 3 * data.getN().size2()) {
431 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
432 "data inconsistency, side %d type %d", side, type);
433 }
434 for (unsigned int gg = 0; gg < data.getN().size1(); ++gg) {
435 for (int dd = 0; dd < 3; dd++) {
436 if ((type == MBTRI && valid_faces3[side]) ||
437 (type == MBEDGE && valid_edges3[side])) {
438 cOords_at_GaussPtF3(gg, dd) += cblas_ddot(
439 nb_dofs / 3, &data.getN(gg)[0], 1, &data.getFieldData()[dd], 3);
440 tAngent1_at_GaussPtF3(gg, dd) +=
441 cblas_ddot(nb_dofs / 3, &data.getDiffN()(gg, 0), 2,
442 &data.getFieldData()[dd], 3);
443 tAngent2_at_GaussPtF3(gg, dd) +=
444 cblas_ddot(nb_dofs / 3, &data.getDiffN()(gg, 1), 2,
445 &data.getFieldData()[dd], 3);
446 } else if ((type == MBTRI && valid_faces4[side]) ||
447 (type == MBEDGE && valid_edges4[side])) {
448 cOords_at_GaussPtF4(gg, dd) += cblas_ddot(
449 nb_dofs / 3, &data.getN(gg)[0], 1, &data.getFieldData()[dd], 3);
450 tAngent1_at_GaussPtF4(gg, dd) +=
451 cblas_ddot(nb_dofs / 3, &data.getDiffN()(gg, 0), 2,
452 &data.getFieldData()[dd], 3);
453 tAngent2_at_GaussPtF4(gg, dd) +=
454 cblas_ddot(nb_dofs / 3, &data.getDiffN()(gg, 1), 2,
455 &data.getFieldData()[dd], 3);
456 }
457 }
458 }
459 } break;
460 default:
461 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
462 }
463
465}
466
469
470 sPin.resize(3, 3);
471 sPin.clear();
472 nOrmals_at_GaussPtF3.resize(tAngent1_at_GaussPtF3.size1(), 3, false);
473 for (unsigned int gg = 0; gg < tAngent1_at_GaussPtF3.size1(); ++gg) {
474 ierr = Spin(&*sPin.data().begin(), &tAngent1_at_GaussPtF3(gg, 0));
475 CHKERRG(ierr);
476 cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1., &*sPin.data().begin(), 3,
477 &tAngent2_at_GaussPtF3(gg, 0), 1, 0.,
478 &nOrmals_at_GaussPtF3(gg, 0), 1);
479 }
480 sPin.clear();
481 nOrmals_at_GaussPtF4.resize(tAngent1_at_GaussPtF4.size1(), 3, false);
482 for (unsigned int gg = 0; gg < tAngent1_at_GaussPtF4.size1(); ++gg) {
483 ierr = Spin(&*sPin.data().begin(), &tAngent1_at_GaussPtF4(gg, 0));
484 CHKERRG(ierr);
485 cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1., &*sPin.data().begin(), 3,
486 &tAngent2_at_GaussPtF4(gg, 0), 1, 0.,
487 &nOrmals_at_GaussPtF4(gg, 0), 1);
488 }
490}
491
493 int side, EntityType type, EntitiesFieldData::EntData &data) {
496
497 if (moab::CN::Dimension(type) != 2)
499
500 if (normalRawPtr == nullptr && normalsAtGaussPtsRawPtr == nullptr)
501 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
502 "Pointer to normal/normals not set");
503
504 bool normal_is_at_gauss_pts = (normalsAtGaussPtsRawPtr != nullptr);
505 if (normal_is_at_gauss_pts)
506 normal_is_at_gauss_pts = (normalsAtGaussPtsRawPtr->size1() != 0);
507
508 auto apply_transform_linear_geometry = [&](auto base, auto nb_gauss_pts,
509 auto nb_base_functions) {
511 const auto &normal = *normalRawPtr;
512 auto t_normal = FTensor::Tensor1<double, 3>{normal[normalShift + 0],
513 normal[normalShift + 1],
514 normal[normalShift + 2]};
515 const auto l02 = t_normal(i) * t_normal(i);
516 auto t_base = data.getFTensor1N<3>(base);
517 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
518 for (int bb = 0; bb != nb_base_functions; ++bb) {
519 const auto v = t_base(0);
520 t_base(i) = (v / l02) * t_normal(i);
521 ++t_base;
522 }
523 }
525 };
526
527 auto apply_transform_nonlinear_geometry = [&](auto base, auto nb_gauss_pts,
528 auto nb_base_functions) {
530 const MatrixDouble &normals_at_pts = *normalsAtGaussPtsRawPtr;
532 &normals_at_pts(0, 0), &normals_at_pts(0, 1), &normals_at_pts(0, 2));
533
534 auto t_base = data.getFTensor1N<3>(base);
535 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
536 const auto l2 = t_normal(i) * t_normal(i);
537 for (int bb = 0; bb != nb_base_functions; ++bb) {
538 const auto v = t_base(0);
539 t_base(i) = (v / l2) * t_normal(i);
540 ++t_base;
541 }
542 ++t_normal;
543 }
545 };
546
547 if (normal_is_at_gauss_pts) {
548 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
549
550 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
551 const auto &base_functions = data.getN(base);
552 const auto nb_gauss_pts = base_functions.size1();
553
554 if (nb_gauss_pts) {
555
556 if (normalsAtGaussPtsRawPtr->size1() != nb_gauss_pts)
557 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
558 "normalsAtGaussPtsRawPtr has inconsistent number of "
559 "integration "
560 "points");
561
562 const auto nb_base_functions = base_functions.size2() / 3;
563 CHKERR apply_transform_nonlinear_geometry(base, nb_gauss_pts,
564 nb_base_functions);
565 }
566 }
567 } else {
568 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
569
570 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
571 const auto &base_functions = data.getN(base);
572 const auto nb_gauss_pts = base_functions.size1();
573
574 if (nb_gauss_pts) {
575 const auto nb_base_functions = base_functions.size2() / 3;
576 CHKERR apply_transform_linear_geometry(base, nb_gauss_pts,
577 nb_base_functions);
578 }
579 }
580 }
581
583}
584
586 int side, EntityType type, EntitiesFieldData::EntData &data) {
588
589 const auto type_dim = moab::CN::Dimension(type);
590 if (type_dim != 1 && type_dim != 2)
592
596
598 &tAngent0[0], &tAngent1[0], &nOrmal[0],
599
600 &tAngent0[1], &tAngent1[1], &nOrmal[1],
601
602 &tAngent0[2], &tAngent1[2], &nOrmal[2]);
603 double det;
606 CHKERR invertTensor3by3(t_m, det, t_inv_m);
607
608 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
609
610 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
611
612 auto &baseN = data.getN(base);
613 auto &diffBaseN = data.getDiffN(base);
614
615 int nb_dofs = baseN.size2() / 3;
616 int nb_gauss_pts = baseN.size1();
617
618 MatrixDouble piola_n(baseN.size1(), baseN.size2());
619 MatrixDouble diff_piola_n(diffBaseN.size1(), diffBaseN.size2());
620
621 if (nb_dofs > 0 && nb_gauss_pts > 0) {
622
624 &baseN(0, HVEC0), &baseN(0, HVEC1), &baseN(0, HVEC2));
626 &diffBaseN(0, HVEC0_0), &diffBaseN(0, HVEC0_1),
627 &diffBaseN(0, HVEC1_0), &diffBaseN(0, HVEC1_1),
628 &diffBaseN(0, HVEC2_0), &diffBaseN(0, HVEC2_1));
629 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_transformed_h_curl(
630 &piola_n(0, HVEC0), &piola_n(0, HVEC1), &piola_n(0, HVEC2));
632 t_transformed_diff_h_curl(
633 &diff_piola_n(0, HVEC0_0), &diff_piola_n(0, HVEC0_1),
634 &diff_piola_n(0, HVEC1_0), &diff_piola_n(0, HVEC1_1),
635 &diff_piola_n(0, HVEC2_0), &diff_piola_n(0, HVEC2_1));
636
637 int cc = 0;
638 if (normalsAtGaussPts.size1() == (unsigned int)nb_gauss_pts) {
639 // HO geometry is set, so jacobian is different at each gauss point
645 &normalsAtGaussPts(0, 2));
646 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
647 CHKERR determinantTensor3by3(t_m_at_pts, det);
648 CHKERR invertTensor3by3(t_m_at_pts, det, t_inv_m);
649 for (int ll = 0; ll != nb_dofs; ll++) {
650 t_transformed_h_curl(i) = t_inv_m(j, i) * t_h_curl(j);
651 t_transformed_diff_h_curl(i, k) =
652 t_inv_m(j, i) * t_diff_h_curl(j, k);
653 ++t_h_curl;
654 ++t_transformed_h_curl;
655 ++t_diff_h_curl;
656 ++t_transformed_diff_h_curl;
657 ++cc;
658 }
659 ++t_m_at_pts;
660 }
661 } else {
662 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
663 for (int ll = 0; ll != nb_dofs; ll++) {
664 t_transformed_h_curl(i) = t_inv_m(j, i) * t_h_curl(j);
665 t_transformed_diff_h_curl(i, k) =
666 t_inv_m(j, i) * t_diff_h_curl(j, k);
667 ++t_h_curl;
668 ++t_transformed_h_curl;
669 ++t_diff_h_curl;
670 ++t_transformed_diff_h_curl;
671 ++cc;
672 }
673 }
674 }
675 if (cc != nb_gauss_pts * nb_dofs)
676 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
677
678 baseN.swap(piola_n);
679 diffBaseN.swap(diff_piola_n);
680 }
681 }
682
684}
685
690
691 int nb_dofs = data.getFieldData().size();
692 if (nb_dofs == 0)
694
695 int nb_gauss_pts = data.getN().size1();
696 tAngent.resize(nb_gauss_pts, 3, false);
697
698 int nb_approx_fun = data.getN().size2();
699 double *diff = &*data.getDiffN().data().begin();
700 double *dofs[] = {&data.getFieldData()[0], &data.getFieldData()[1],
701 &data.getFieldData()[2]};
702
703 tAngent.resize(nb_gauss_pts, 3, false);
704
705 switch (type) {
706 case MBVERTEX:
707 for (int dd = 0; dd != 3; dd++) {
708 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
709 tAngent(gg, dd) = cblas_ddot(2, diff, 1, dofs[dd], 3);
710 }
711 }
712 break;
713 case MBEDGE:
714 if (nb_dofs % 3) {
715 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
716 "Approximated field should be rank 3, i.e. vector in 3d space");
717 }
718 for (int dd = 0; dd != 3; dd++) {
719 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
720 tAngent(gg, dd) +=
721 cblas_ddot(nb_dofs / 3, &diff[gg * nb_approx_fun], 1, dofs[dd], 3);
722 }
723 }
724 break;
725 default:
726 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
727 "This operator can calculate tangent vector only on edge");
728 }
729
731}
732
734 int side, EntityType type, EntitiesFieldData::EntData &data) {
736
737 if (type != MBEDGE)
739
742 &tAngent[0], &tAngent[1], &tAngent[2]);
743 const double l0 = t_m(i) * t_m(i);
744
745 auto get_base_at_pts = [&](auto base) {
747 &data.getN(base)(0, HVEC0), &data.getN(base)(0, HVEC1),
748 &data.getN(base)(0, HVEC2));
749 return t_h_curl;
750 };
751
752 auto get_tangent_at_pts = [&]() {
754 &tangentAtGaussPt(0, 0), &tangentAtGaussPt(0, 1),
755 &tangentAtGaussPt(0, 2));
756 return t_m_at_pts;
757 };
758
759 auto calculate_squared_edge_length = [&]() {
760 std::vector<double> l1;
761 int nb_gauss_pts = tangentAtGaussPt.size1();
762 if (nb_gauss_pts) {
763 l1.resize(nb_gauss_pts);
764 auto t_m_at_pts = get_tangent_at_pts();
765 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
766 l1[gg] = t_m_at_pts(i) * t_m_at_pts(i);
767 ++t_m_at_pts;
768 }
769 }
770 return l1;
771 };
772
773 auto l1 = calculate_squared_edge_length();
774
775 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
776
777 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
778 const size_t nb_gauss_pts = data.getN(base).size1();
779 const size_t nb_dofs = data.getN(base).size2() / 3;
780 if (nb_gauss_pts && nb_dofs) {
781 auto t_h_curl = get_base_at_pts(base);
782 int cc = 0;
783 if (tangentAtGaussPt.size1() == nb_gauss_pts) {
784 auto t_m_at_pts = get_tangent_at_pts();
785 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
786 const double l0 = l1[gg];
787 for (int ll = 0; ll != nb_dofs; ll++) {
788 const double val = t_h_curl(0);
789 const double a = val / l0;
790 t_h_curl(i) = t_m_at_pts(i) * a;
791 ++t_h_curl;
792 ++cc;
793 }
794 ++t_m_at_pts;
795 }
796 } else {
797 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
798 for (int ll = 0; ll != nb_dofs; ll++) {
799 const double val = t_h_curl(0);
800 const double a = val / l0;
801 t_h_curl(i) = t_m(i) * a;
802 ++t_h_curl;
803 ++cc;
804 }
805 }
806 }
807
808 if (cc != nb_gauss_pts * nb_dofs)
809 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
810 }
811 }
812
814}
815
816template <>
817template <>
820 double *ptr = &*data.data().begin();
821 return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2], 3);
822}
823
824template <>
825template <>
827OpGetDataAndGradient<3, 3>::getGradAtGaussPtsTensor<3, 3>(MatrixDouble &data) {
828 double *ptr = &*data.data().begin();
829 return FTensor::Tensor2<double *, 3, 3>(ptr, &ptr[1], &ptr[2], &ptr[3],
830 &ptr[4], &ptr[5], &ptr[6], &ptr[7],
831 &ptr[8], 9);
832}
833
834template <>
836 int side, EntityType type, EntitiesFieldData::EntData &data) {
838 if (data.getBase() == NOBASE)
840 const unsigned int nb_gauss_pts = data.getN().size1();
841 const unsigned int nb_base_functions = data.getN().size2();
842 const unsigned int nb_dofs = data.getFieldData().size();
843 if (!nb_dofs)
845 auto t_n = data.getFTensor0N();
846 auto t_val = getValAtGaussPtsTensor<3>(dataAtGaussPts);
847 auto t_grad = getGradAtGaussPtsTensor<3, 3>(dataGradAtGaussPts);
850 if (type == MBVERTEX &&
851 data.getDiffN().data().size() == 3 * nb_base_functions) {
852 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
853 auto t_data = data.getFTensor1FieldData<3>();
854 auto t_diff_n = data.getFTensor1DiffN<3>();
855 unsigned int bb = 0;
856 for (; bb != nb_dofs / 3; ++bb) {
857 t_val(i) += t_data(i) * t_n;
858 t_grad(i, j) += t_data(i) * t_diff_n(j);
859 ++t_n;
860 ++t_diff_n;
861 ++t_data;
862 }
863 ++t_val;
864 ++t_grad;
865 for (; bb != nb_base_functions; ++bb) {
866 ++t_n;
867 }
868 }
869 } else {
870 auto t_diff_n = data.getFTensor1DiffN<3>();
871 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
872 auto t_data = data.getFTensor1FieldData<3>();
873 unsigned int bb = 0;
874 for (; bb != nb_dofs / 3; ++bb) {
875 t_val(i) += t_data(i) * t_n;
876 t_grad(i, j) += t_data(i) * t_diff_n(j);
877 ++t_n;
878 ++t_diff_n;
879 ++t_data;
880 }
881 ++t_val;
882 ++t_grad;
883 for (; bb != nb_base_functions; ++bb) {
884 ++t_n;
885 ++t_diff_n;
886 }
887 }
888 }
890}
891
892template <>
894 int side, EntityType type, EntitiesFieldData::EntData &data) {
896 const unsigned int nb_gauss_pts = data.getN().size1();
897 const unsigned int nb_base_functions = data.getN().size2();
898 // bool constant_diff = false;
899 const unsigned int nb_dofs = data.getFieldData().size();
900 auto t_n = data.getFTensor0N();
902 FTensor::Tensor0<double *>(&*dataAtGaussPts.data().begin(), 1);
903 double *ptr = &*dataGradAtGaussPts.data().begin();
905 &ptr[2]);
907 if (type == MBVERTEX &&
908 data.getDiffN().data().size() == 3 * nb_base_functions) {
909 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
910 auto t_data = data.getFTensor0FieldData();
911 auto t_diff_n = data.getFTensor1DiffN<3>();
912 unsigned int bb = 0;
913 for (; bb != nb_dofs / 3; ++bb) {
914 t_val += t_data * t_n;
915 t_grad(i) += t_data * t_diff_n(i);
916 ++t_n;
917 ++t_diff_n;
918 ++t_data;
919 }
920 ++t_val;
921 ++t_grad;
922 for (; bb != nb_base_functions; ++bb) {
923 ++t_n;
924 }
925 }
926 } else {
927 auto t_diff_n = data.getFTensor1DiffN<3>();
928 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
929 auto t_data = data.getFTensor0FieldData();
930 unsigned int bb = 0;
931 for (; bb != nb_dofs / 3; ++bb) {
932 t_val = t_data * t_n;
933 t_grad(i) += t_data * t_diff_n(i);
934 ++t_n;
935 ++t_diff_n;
936 ++t_data;
937 }
938 ++t_val;
939 ++t_grad;
940 for (; bb != nb_base_functions; ++bb) {
941 ++t_n;
942 ++t_diff_n;
943 }
944 }
945 }
947}
948} // namespace MoFEM
constexpr double a
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ NOBASE
Definition: definitions.h:59
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ HVEC0
Definition: definitions.h:186
@ HVEC1
Definition: definitions.h:186
@ HVEC2
Definition: definitions.h:186
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
@ HVEC1_1
Definition: definitions.h:196
@ HVEC0_1
Definition: definitions.h:195
@ HVEC1_0
Definition: definitions.h:193
@ HVEC2_1
Definition: definitions.h:197
@ HVEC1_2
Definition: definitions.h:199
@ HVEC2_2
Definition: definitions.h:200
@ HVEC2_0
Definition: definitions.h:194
@ HVEC0_2
Definition: definitions.h:198
@ HVEC0_0
Definition: definitions.h:192
#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
PetscErrorCode Spin(double *spinOmega, double *vecOmega)
calculate spin matrix from vector
Definition: fem_tools.c:546
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1608
MoFEMErrorCode invertTensor3by3< 3, double, ublas::row_major, DoubleAllocator >(MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data)
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1557
constexpr IntegrationType I
constexpr AssemblyType A
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
virtual MoFEMErrorCode opLhs(EntitiesFieldData &row_data, EntitiesFieldData &col_data)
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
virtual MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
bool getSymm() const
Get if operator uses symmetry of DOFs or not.
DataOperator(const bool symm=true)
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBDiffNByOrderArray()
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
FieldApproximationBase & getBase()
Get approximation base.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBDiffNMap()
get hash map of direvarives base function for BB base, key is a field name
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
data structure for finite element entity
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Get field values and gradients at Gauss points.
MoFEMErrorCode calculateValAndGrad(int side, EntityType type, EntitiesFieldData::EntData &data)
Calculate gradient and values at integration points.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Tensor2< double *, 3, 3 > tJac
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
int normalShift
Shift in vector for linear geometry.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Tensor2< double *, 3, 3 > tInvJac
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'j', 3 > j
MatrixDouble diffNinvJac
FTensor::Index< 'i', 3 > i
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Tensor2< double *, 3, 3 > tInvJac
FTensor::Index< 'j', 3 > j
FTensor::Tensor2< double *, 3, 3 > tInvJac
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k