v0.15.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
169MoFEMErrorCode OpSetInvJacH1::doWork(int side, EntityType type,
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
222OpSetInvJacHdivAndHcurl::doWork(int side, EntityType type,
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 *inv_diff_n_ptr = &*diffHdivInvJac.data().begin();
243 inv_diff_n_ptr, &inv_diff_n_ptr[HVEC0_1], &inv_diff_n_ptr[HVEC0_2],
244
245 &inv_diff_n_ptr[HVEC1_0], &inv_diff_n_ptr[HVEC1_1],
246 &inv_diff_n_ptr[HVEC1_2],
247
248 &inv_diff_n_ptr[HVEC2_0], &inv_diff_n_ptr[HVEC2_1],
249 &inv_diff_n_ptr[HVEC2_2]);
250
251 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
252 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
253 t_inv_diff_n(k, i) = t_diff_n(k, j) * tInvJac(j, i);
254 ++t_diff_n;
255 ++t_inv_diff_n;
256 }
257 }
258
259 data.getDiffN(base).swap(diffHdivInvJac);
260 }
261
263}
264
266 int side, EntityType type, EntitiesFieldData::EntData &data) {
268
269 if (CN::Dimension(type) > 1) {
270
271 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
272
273 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
274
275 const unsigned int nb_base_functions = data.getN(base).size2() / 3;
276 if (!nb_base_functions)
277 continue;
278
279 const unsigned int nb_gauss_pts = data.getN(base).size1();
280 double const a = 1. / vOlume;
281
282 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
283 if (data.getN(base).size2() > 0) {
284 auto t_n = data.getFTensor1N<3>(base);
285 double *t_transformed_n_ptr = &*piolaN.data().begin();
287 t_transformed_n_ptr, // HVEC0
288 &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
289 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
290 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
291 t_transformed_n(i) = a * (tJac(i, k) * t_n(k));
292 ++t_n;
293 ++t_transformed_n;
294 }
295 }
296 data.getN(base).swap(piolaN);
297 }
298
299 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
300 if (data.getDiffN(base).size2() > 0) {
301 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
302 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
304 t_transformed_diff_n(t_transformed_diff_n_ptr,
305 &t_transformed_diff_n_ptr[HVEC0_1],
306 &t_transformed_diff_n_ptr[HVEC0_2],
307 &t_transformed_diff_n_ptr[HVEC1_0],
308 &t_transformed_diff_n_ptr[HVEC1_1],
309 &t_transformed_diff_n_ptr[HVEC1_2],
310 &t_transformed_diff_n_ptr[HVEC2_0],
311 &t_transformed_diff_n_ptr[HVEC2_1],
312 &t_transformed_diff_n_ptr[HVEC2_2]);
313 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
314 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
315 t_transformed_diff_n(i, k) = a * tJac(i, j) * t_diff_n(j, k);
316 ++t_diff_n;
317 ++t_transformed_diff_n;
318 }
319 }
320 data.getDiffN(base).swap(piolaDiffN);
321 }
322 }
323 }
324
326}
327
329OpSetCovariantPiolaTransform::doWork(int side, EntityType type,
332
333 if (type == MBVERTEX)
335
336 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
337
338 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
339
340 const unsigned int nb_base_functions = data.getN(base).size2() / 3;
341 if (!nb_base_functions)
342 continue;
343
344 const unsigned int nb_gauss_pts = data.getN(base).size1();
345 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
346 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
347
348 auto t_n = data.getFTensor1N<3>(base);
349 double *t_transformed_n_ptr = &*piolaN.data().begin();
351 t_transformed_n_ptr, &t_transformed_n_ptr[HVEC1],
352 &t_transformed_n_ptr[HVEC2]);
353 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
354 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
355 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_transformed_diff_n(
356 t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[HVEC0_1],
357 &t_transformed_diff_n_ptr[HVEC0_2], &t_transformed_diff_n_ptr[HVEC1_0],
358 &t_transformed_diff_n_ptr[HVEC1_1], &t_transformed_diff_n_ptr[HVEC1_2],
359 &t_transformed_diff_n_ptr[HVEC2_0], &t_transformed_diff_n_ptr[HVEC2_1],
360 &t_transformed_diff_n_ptr[HVEC2_2]);
361
362 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
363 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
364 t_transformed_n(i) = tInvJac(k, i) * t_n(k);
365 ++t_n;
366 ++t_transformed_n;
367 }
368 }
369
370 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
371 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
372 t_transformed_diff_n(i, k) = tInvJac(j, i) * t_diff_n(j, k);
373 ++t_diff_n;
374 ++t_transformed_diff_n;
375 }
376 }
377
378 data.getN(base).swap(piolaN);
379 data.getDiffN(base).swap(piolaDiffN);
380 }
381
382 // data.getBase() = base;
383
385}
386
388OpGetCoordsAndNormalsOnPrism::doWork(int side, EntityType type,
391
392 if (data.getFieldData().size() == 0)
394 const int valid_edges3[] = {1, 1, 1, 0, 0, 0, 0, 0, 0};
395 const int valid_faces3[] = {0, 0, 0, 1, 0, 0, 0, 0, 0};
396 const int valid_edges4[] = {0, 0, 0, 0, 0, 0, 1, 1, 1};
397 const int valid_faces4[] = {0, 0, 0, 0, 1, 0, 0, 0, 0};
398
399 if (type == MBEDGE) {
400 if (!valid_edges3[side] || valid_edges4[side])
402 } else if (type == MBTRI) {
403 if (!valid_faces3[side] || valid_faces4[side])
405 }
406
407 switch (type) {
408 case MBVERTEX: {
409 for (unsigned int gg = 0; gg < data.getN().size1(); ++gg) {
410 for (int dd = 0; dd < 3; dd++) {
411 cOords_at_GaussPtF3(gg, dd) =
412 cblas_ddot(3, &data.getN(gg)[0], 1, &data.getFieldData()[dd], 3);
413 tAngent1_at_GaussPtF3(gg, dd) = cblas_ddot(
414 3, &data.getDiffN()(gg, 0), 2, &data.getFieldData()[dd], 3);
415 tAngent2_at_GaussPtF3(gg, dd) = cblas_ddot(
416 3, &data.getDiffN()(gg, 1), 2, &data.getFieldData()[dd], 3);
417 cOords_at_GaussPtF4(gg, dd) = cblas_ddot(
418 3, &data.getN(gg)[0], 1, &data.getFieldData()[9 + dd], 3);
419 tAngent1_at_GaussPtF4(gg, dd) = cblas_ddot(
420 3, &data.getDiffN()(gg, 6 + 0), 2, &data.getFieldData()[9 + dd], 3);
421 tAngent2_at_GaussPtF4(gg, dd) = cblas_ddot(
422 3, &data.getDiffN()(gg, 6 + 1), 2, &data.getFieldData()[9 + dd], 3);
423 }
424 }
425 } break;
426 case MBEDGE:
427 case MBTRI: {
428 if (2 * data.getN().size2() != data.getDiffN().size2()) {
429 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
430 }
431 unsigned int nb_dofs = data.getFieldData().size();
432 if (nb_dofs % 3 != 0) {
433 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
434 }
435 if (nb_dofs > 3 * data.getN().size2()) {
436 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
437 "data inconsistency, side %d type %d", side, type);
438 }
439 for (unsigned int gg = 0; gg < data.getN().size1(); ++gg) {
440 for (int dd = 0; dd < 3; dd++) {
441 if ((type == MBTRI && valid_faces3[side]) ||
442 (type == MBEDGE && valid_edges3[side])) {
443 cOords_at_GaussPtF3(gg, dd) += cblas_ddot(
444 nb_dofs / 3, &data.getN(gg)[0], 1, &data.getFieldData()[dd], 3);
445 tAngent1_at_GaussPtF3(gg, dd) +=
446 cblas_ddot(nb_dofs / 3, &data.getDiffN()(gg, 0), 2,
447 &data.getFieldData()[dd], 3);
448 tAngent2_at_GaussPtF3(gg, dd) +=
449 cblas_ddot(nb_dofs / 3, &data.getDiffN()(gg, 1), 2,
450 &data.getFieldData()[dd], 3);
451 } else if ((type == MBTRI && valid_faces4[side]) ||
452 (type == MBEDGE && valid_edges4[side])) {
453 cOords_at_GaussPtF4(gg, dd) += cblas_ddot(
454 nb_dofs / 3, &data.getN(gg)[0], 1, &data.getFieldData()[dd], 3);
455 tAngent1_at_GaussPtF4(gg, dd) +=
456 cblas_ddot(nb_dofs / 3, &data.getDiffN()(gg, 0), 2,
457 &data.getFieldData()[dd], 3);
458 tAngent2_at_GaussPtF4(gg, dd) +=
459 cblas_ddot(nb_dofs / 3, &data.getDiffN()(gg, 1), 2,
460 &data.getFieldData()[dd], 3);
461 }
462 }
463 }
464 } break;
465 default:
466 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
467 }
468
470}
471
474
475 sPin.resize(3, 3);
476 sPin.clear();
477 nOrmals_at_GaussPtF3.resize(tAngent1_at_GaussPtF3.size1(), 3, false);
478 for (unsigned int gg = 0; gg < tAngent1_at_GaussPtF3.size1(); ++gg) {
479 ierr = Spin(&*sPin.data().begin(), &tAngent1_at_GaussPtF3(gg, 0));
480 CHKERRG(ierr);
481 cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1., &*sPin.data().begin(), 3,
482 &tAngent2_at_GaussPtF3(gg, 0), 1, 0.,
483 &nOrmals_at_GaussPtF3(gg, 0), 1);
484 }
485 sPin.clear();
486 nOrmals_at_GaussPtF4.resize(tAngent1_at_GaussPtF4.size1(), 3, false);
487 for (unsigned int gg = 0; gg < tAngent1_at_GaussPtF4.size1(); ++gg) {
488 ierr = Spin(&*sPin.data().begin(), &tAngent1_at_GaussPtF4(gg, 0));
489 CHKERRG(ierr);
490 cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1., &*sPin.data().begin(), 3,
491 &tAngent2_at_GaussPtF4(gg, 0), 1, 0.,
492 &nOrmals_at_GaussPtF4(gg, 0), 1);
493 }
495}
496
498 int side, EntityType type, EntitiesFieldData::EntData &data) {
499 FTensor::Index<'i', 3> i;
501
502 if (moab::CN::Dimension(type) != 2)
504
505 if (normalRawPtr == nullptr && normalsAtGaussPtsRawPtr == nullptr)
506 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
507 "Pointer to normal/normals not set");
508
509 bool normal_is_at_gauss_pts = (normalsAtGaussPtsRawPtr != nullptr);
510 if (normal_is_at_gauss_pts)
511 normal_is_at_gauss_pts = (normalsAtGaussPtsRawPtr->size1() != 0);
512
513 auto apply_transform_linear_geometry = [&](auto base, auto nb_gauss_pts,
514 auto nb_base_functions) {
516 const auto &normal = *normalRawPtr;
517 auto t_normal = FTensor::Tensor1<double, 3>{normal[normalShift + 0],
518 normal[normalShift + 1],
519 normal[normalShift + 2]};
520 const auto l02 = t_normal(i) * t_normal(i);
521 auto t_base = data.getFTensor1N<3>(base);
522 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
523 for (int bb = 0; bb != nb_base_functions; ++bb) {
524 const auto v = t_base(0);
525 t_base(i) = (v / l02) * t_normal(i);
526 ++t_base;
527 }
528 }
530 };
531
532 auto apply_transform_nonlinear_geometry = [&](auto base, auto nb_gauss_pts,
533 auto nb_base_functions) {
535 const MatrixDouble &normals_at_pts = *normalsAtGaussPtsRawPtr;
537 &normals_at_pts(0, 0), &normals_at_pts(0, 1), &normals_at_pts(0, 2));
538
539 auto t_base = data.getFTensor1N<3>(base);
540 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
541 const auto l2 = t_normal(i) * t_normal(i);
542 for (int bb = 0; bb != nb_base_functions; ++bb) {
543 const auto v = t_base(0);
544 t_base(i) = (v / l2) * t_normal(i);
545 ++t_base;
546 }
547 ++t_normal;
548 }
550 };
551
552 if (normal_is_at_gauss_pts) {
553 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
554
555 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
556 const auto &base_functions = data.getN(base);
557 const auto nb_gauss_pts = base_functions.size1();
558
559 if (nb_gauss_pts) {
560
561 if (normalsAtGaussPtsRawPtr->size1() != nb_gauss_pts)
562 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
563 "normalsAtGaussPtsRawPtr has inconsistent number of "
564 "integration "
565 "points");
566
567 const auto nb_base_functions = base_functions.size2() / 3;
568 CHKERR apply_transform_nonlinear_geometry(base, nb_gauss_pts,
569 nb_base_functions);
570 }
571 }
572 } else {
573 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
574
575 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
576 const auto &base_functions = data.getN(base);
577 const auto nb_gauss_pts = base_functions.size1();
578
579 if (nb_gauss_pts) {
580 const auto nb_base_functions = base_functions.size2() / 3;
581 CHKERR apply_transform_linear_geometry(base, nb_gauss_pts,
582 nb_base_functions);
583 }
584 }
585 }
586
588}
589
591 int side, EntityType type, EntitiesFieldData::EntData &data) {
593
594 const auto type_dim = moab::CN::Dimension(type);
595 if (type_dim != 1 && type_dim != 2)
597
598 FTensor::Index<'i', 3> i;
599 FTensor::Index<'j', 3> j;
600 FTensor::Index<'k', 2> k;
601
603 &tAngent0[0], &tAngent1[0], &nOrmal[0],
604
605 &tAngent0[1], &tAngent1[1], &nOrmal[1],
606
607 &tAngent0[2], &tAngent1[2], &nOrmal[2]);
608 double det;
611 CHKERR invertTensor3by3(t_m, det, t_inv_m);
612
613 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
614
615 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
616
617 auto &baseN = data.getN(base);
618 auto &diffBaseN = data.getDiffN(base);
619
620 int nb_dofs = baseN.size2() / 3;
621 int nb_gauss_pts = baseN.size1();
622
623 MatrixDouble piola_n(baseN.size1(), baseN.size2());
624 MatrixDouble diff_piola_n(diffBaseN.size1(), diffBaseN.size2());
625
626 if (nb_dofs > 0 && nb_gauss_pts > 0) {
627
629 &baseN(0, HVEC0), &baseN(0, HVEC1), &baseN(0, HVEC2));
631 &diffBaseN(0, HVEC0_0), &diffBaseN(0, HVEC0_1),
632 &diffBaseN(0, HVEC1_0), &diffBaseN(0, HVEC1_1),
633 &diffBaseN(0, HVEC2_0), &diffBaseN(0, HVEC2_1));
634 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_transformed_h_curl(
635 &piola_n(0, HVEC0), &piola_n(0, HVEC1), &piola_n(0, HVEC2));
637 t_transformed_diff_h_curl(
638 &diff_piola_n(0, HVEC0_0), &diff_piola_n(0, HVEC0_1),
639 &diff_piola_n(0, HVEC1_0), &diff_piola_n(0, HVEC1_1),
640 &diff_piola_n(0, HVEC2_0), &diff_piola_n(0, HVEC2_1));
641
642 int cc = 0;
643 if (normalsAtGaussPts.size1() == (unsigned int)nb_gauss_pts) {
644 // HO geometry is set, so jacobian is different at each gauss point
650 &normalsAtGaussPts(0, 2));
651 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
652 CHKERR determinantTensor3by3(t_m_at_pts, det);
653 CHKERR invertTensor3by3(t_m_at_pts, det, t_inv_m);
654 for (int ll = 0; ll != nb_dofs; ll++) {
655 t_transformed_h_curl(i) = t_inv_m(j, i) * t_h_curl(j);
656 t_transformed_diff_h_curl(i, k) =
657 t_inv_m(j, i) * t_diff_h_curl(j, k);
658 ++t_h_curl;
659 ++t_transformed_h_curl;
660 ++t_diff_h_curl;
661 ++t_transformed_diff_h_curl;
662 ++cc;
663 }
664 ++t_m_at_pts;
665 }
666 } else {
667 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
668 for (int ll = 0; ll != nb_dofs; ll++) {
669 t_transformed_h_curl(i) = t_inv_m(j, i) * t_h_curl(j);
670 t_transformed_diff_h_curl(i, k) =
671 t_inv_m(j, i) * t_diff_h_curl(j, k);
672 ++t_h_curl;
673 ++t_transformed_h_curl;
674 ++t_diff_h_curl;
675 ++t_transformed_diff_h_curl;
676 ++cc;
677 }
678 }
679 }
680 if (cc != nb_gauss_pts * nb_dofs)
681 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
682
683 baseN.swap(piola_n);
684 diffBaseN.swap(diff_piola_n);
685 }
686 }
687
689}
690
692OpGetHOTangentOnEdge::doWork(int side, EntityType type,
695
696 int nb_dofs = data.getFieldData().size();
697 if (nb_dofs == 0)
699
700 int nb_gauss_pts = data.getN().size1();
701 tAngent.resize(nb_gauss_pts, 3, false);
702
703 int nb_approx_fun = data.getN().size2();
704 double *diff = &*data.getDiffN().data().begin();
705 double *dofs[] = {&data.getFieldData()[0], &data.getFieldData()[1],
706 &data.getFieldData()[2]};
707
708 tAngent.resize(nb_gauss_pts, 3, false);
709
710 switch (type) {
711 case MBVERTEX:
712 for (int dd = 0; dd != 3; dd++) {
713 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
714 tAngent(gg, dd) = cblas_ddot(2, diff, 1, dofs[dd], 3);
715 }
716 }
717 break;
718 case MBEDGE:
719 if (nb_dofs % 3) {
720 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
721 "Approximated field should be rank 3, i.e. vector in 3d space");
722 }
723 for (int dd = 0; dd != 3; dd++) {
724 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
725 tAngent(gg, dd) +=
726 cblas_ddot(nb_dofs / 3, &diff[gg * nb_approx_fun], 1, dofs[dd], 3);
727 }
728 }
729 break;
730 default:
731 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
732 "This operator can calculate tangent vector only on edge");
733 }
734
736}
737
739 int side, EntityType type, EntitiesFieldData::EntData &data) {
741
742 if (type != MBEDGE)
744
745 FTensor::Index<'i', 3> i;
747 &tAngent[0], &tAngent[1], &tAngent[2]);
748 const double l0 = t_m(i) * t_m(i);
749
750 auto get_base_at_pts = [&](auto base) {
752 &data.getN(base)(0, HVEC0), &data.getN(base)(0, HVEC1),
753 &data.getN(base)(0, HVEC2));
754 return t_h_curl;
755 };
756
757 auto get_tangent_at_pts = [&]() {
759 &tangentAtGaussPt(0, 0), &tangentAtGaussPt(0, 1),
760 &tangentAtGaussPt(0, 2));
761 return t_m_at_pts;
762 };
763
764 auto calculate_squared_edge_length = [&]() {
765 std::vector<double> l1;
766 int nb_gauss_pts = tangentAtGaussPt.size1();
767 if (nb_gauss_pts) {
768 l1.resize(nb_gauss_pts);
769 auto t_m_at_pts = get_tangent_at_pts();
770 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
771 l1[gg] = t_m_at_pts(i) * t_m_at_pts(i);
772 ++t_m_at_pts;
773 }
774 }
775 return l1;
776 };
777
778 auto l1 = calculate_squared_edge_length();
779
780 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
781
782 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
783 const size_t nb_gauss_pts = data.getN(base).size1();
784 const size_t nb_dofs = data.getN(base).size2() / 3;
785 if (nb_gauss_pts && nb_dofs) {
786 auto t_h_curl = get_base_at_pts(base);
787 int cc = 0;
788 if (tangentAtGaussPt.size1() == nb_gauss_pts) {
789 auto t_m_at_pts = get_tangent_at_pts();
790 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
791 const double l0 = l1[gg];
792 for (int ll = 0; ll != nb_dofs; ll++) {
793 const double val = t_h_curl(0);
794 const double a = val / l0;
795 t_h_curl(i) = t_m_at_pts(i) * a;
796 ++t_h_curl;
797 ++cc;
798 }
799 ++t_m_at_pts;
800 }
801 } else {
802 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
803 for (int ll = 0; ll != nb_dofs; ll++) {
804 const double val = t_h_curl(0);
805 const double a = val / l0;
806 t_h_curl(i) = t_m(i) * a;
807 ++t_h_curl;
808 ++cc;
809 }
810 }
811 }
812
813 if (cc != nb_gauss_pts * nb_dofs)
814 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
815 }
816 }
817
819}
820
821template <>
822template <>
825 double *ptr = &*data.data().begin();
826 return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2], 3);
827}
828
829template <>
830template <>
833 double *ptr = &*data.data().begin();
834 return FTensor::Tensor2<double *, 3, 3>(ptr, &ptr[1], &ptr[2], &ptr[3],
835 &ptr[4], &ptr[5], &ptr[6], &ptr[7],
836 &ptr[8], 9);
837}
838
839template <>
841 int side, EntityType type, EntitiesFieldData::EntData &data) {
843 if (data.getBase() == NOBASE)
845 const unsigned int nb_gauss_pts = data.getN().size1();
846 const unsigned int nb_base_functions = data.getN().size2();
847 const unsigned int nb_dofs = data.getFieldData().size();
848 if (!nb_dofs)
850 auto t_n = data.getFTensor0N();
851 auto t_val = getValAtGaussPtsTensor<3>(dataAtGaussPts);
852 auto t_grad = getGradAtGaussPtsTensor<3, 3>(dataGradAtGaussPts);
853 FTensor::Index<'i', 3> i;
854 FTensor::Index<'j', 3> j;
855 if (type == MBVERTEX &&
856 data.getDiffN().data().size() == 3 * nb_base_functions) {
857 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
858 auto t_data = data.getFTensor1FieldData<3>();
859 auto t_diff_n = data.getFTensor1DiffN<3>();
860 unsigned int bb = 0;
861 for (; bb != nb_dofs / 3; ++bb) {
862 t_val(i) += t_data(i) * t_n;
863 t_grad(i, j) += t_data(i) * t_diff_n(j);
864 ++t_n;
865 ++t_diff_n;
866 ++t_data;
867 }
868 ++t_val;
869 ++t_grad;
870 for (; bb != nb_base_functions; ++bb) {
871 ++t_n;
872 }
873 }
874 } else {
875 auto t_diff_n = data.getFTensor1DiffN<3>();
876 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
877 auto t_data = data.getFTensor1FieldData<3>();
878 unsigned int bb = 0;
879 for (; bb != nb_dofs / 3; ++bb) {
880 t_val(i) += t_data(i) * t_n;
881 t_grad(i, j) += t_data(i) * t_diff_n(j);
882 ++t_n;
883 ++t_diff_n;
884 ++t_data;
885 }
886 ++t_val;
887 ++t_grad;
888 for (; bb != nb_base_functions; ++bb) {
889 ++t_n;
890 ++t_diff_n;
891 }
892 }
893 }
895}
896
897template <>
899 int side, EntityType type, EntitiesFieldData::EntData &data) {
901 const unsigned int nb_gauss_pts = data.getN().size1();
902 const unsigned int nb_base_functions = data.getN().size2();
903 // bool constant_diff = false;
904 const unsigned int nb_dofs = data.getFieldData().size();
905 auto t_n = data.getFTensor0N();
907 FTensor::Tensor0<double *>(&*dataAtGaussPts.data().begin(), 1);
908 double *ptr = &*dataGradAtGaussPts.data().begin();
910 &ptr[2]);
911 FTensor::Index<'i', 3> i;
912 if (type == MBVERTEX &&
913 data.getDiffN().data().size() == 3 * nb_base_functions) {
914 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
915 auto t_data = data.getFTensor0FieldData();
916 auto t_diff_n = data.getFTensor1DiffN<3>();
917 unsigned int bb = 0;
918 for (; bb != nb_dofs / 3; ++bb) {
919 t_val += t_data * t_n;
920 t_grad(i) += t_data * t_diff_n(i);
921 ++t_n;
922 ++t_diff_n;
923 ++t_data;
924 }
925 ++t_val;
926 ++t_grad;
927 for (; bb != nb_base_functions; ++bb) {
928 ++t_n;
929 }
930 }
931 } else {
932 auto t_diff_n = data.getFTensor1DiffN<3>();
933 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
934 auto t_data = data.getFTensor0FieldData();
935 unsigned int bb = 0;
936 for (; bb != nb_dofs / 3; ++bb) {
937 t_val = t_data * t_n;
938 t_grad(i) += t_data * t_diff_n(i);
939 ++t_n;
940 ++t_diff_n;
941 ++t_data;
942 }
943 ++t_val;
944 ++t_grad;
945 for (; bb != nb_base_functions; ++bb) {
946 ++t_n;
947 ++t_diff_n;
948 }
949 }
950 }
952}
953} // namespace MoFEM
constexpr double a
T data[Tensor_Dim]
FieldApproximationBase
approximation base
Definition definitions.h:58
@ LASTBASE
Definition definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
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()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ HVEC0
@ HVEC1
@ HVEC2
@ 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()
@ HVEC1_1
@ HVEC0_1
@ HVEC1_0
@ HVEC2_1
@ HVEC1_2
@ HVEC2_2
@ HVEC2_0
@ HVEC0_2
@ HVEC0_0
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PetscErrorCode Spin(double *spinOmega, double *vecOmega)
calculate spin matrix from vector
Definition fem_tools.c:546
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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.
MoFEMErrorCode invertTensor3by3< 3, double, ublas::row_major, DoubleAllocator >(MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data)
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
constexpr IntegrationType I
constexpr AssemblyType A
FTensor::Index< 'm', 3 > m
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 field data coefficients.
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBDiffNMap()
get hash map of derivatives 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.
FTensor::Tensor2< double *, R, D > getGradAtGaussPtsTensor(MatrixDouble &data)
MoFEMErrorCode calculateValAndGrad(int side, EntityType type, EntitiesFieldData::EntData &data)
Calculate gradient and values at integration points.
FTensor::Tensor1< double *, R > getValAtGaussPtsTensor(MatrixDouble &data)
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.
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 > tJac
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.
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::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