13 FieldSpace space, boost::shared_ptr<MatrixDouble> inv_jac_ptr)
29 "This operator can be used only with element which faces");
32 return applyTransform<2, 2, 2, 2>(diff_n);
35 if (!(type == MBVERTEX && sPace ==
L2)) {
45 CHKERR apply_transform(*(
m.second));
50 CHKERR apply_transform(*ptr);
64 "This operator can be used only with element which face");
67 return applyTransform<2, 3, 3, 3>(diff_n);
70 if (!(type == MBVERTEX && sPace ==
L2)) {
79 CHKERR apply_transform(*(
m.second));
84 CHKERR apply_transform(*ptr);
98 "This operator can be used only with element which faces");
102 size_t nb_functions = diff2_n.size2() / 4;
104 size_t nb_gauss_pts = diff2_n.size1();
107 if (nb_gauss_pts != getGaussPts().size2())
109 "Wrong number of Gauss Pts");
112 diffNinvJac.resize(diff2_n.size1(), diff2_n.size2(),
false);
120 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
121 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
122 for (
size_t dd = 0; dd != nb_functions; ++dd) {
124 t_inv_jac(
k,
i) * (t_inv_jac(
l,
j) * t_diff2_n_ref(
k,
l));
130 diff2_n.swap(diffNinvJac);
135 if (!(type == MBVERTEX && sPace ==
L2)) {
140 data.
getN(base, BaseDerivatives::SecondDerivative));
143 auto error = [&](
auto &
m) {
147 "Operator do not work for Bernstein-Bezier. This "
157 CHKERR apply_transform(*(
m.second));
164 CHKERR apply_transform(*ptr);
177 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
182 "This operator can be used only with element which is triangle");
192 const unsigned int nb_base_functions = data.
getDiffN(base).size2() / 6;
193 if (nb_base_functions) {
194 const unsigned int nb_gauss_pts = data.
getDiffN(base).size1();
196 diffHcurlInvJac.resize(nb_gauss_pts, data.
getDiffN(base).size2(),
false);
199 double *t_inv_diff_n_ptr = &*diffHcurlInvJac.data().begin();
201 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[
HVEC0_1],
207 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
208 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
209 for (
unsigned int bb = 0; bb != nb_base_functions; bb++) {
210 t_inv_diff_n(
i,
j) = t_diff_n(
i,
k) * t_inv_jac(
k,
j);
216 data.
getDiffN(base).swap(diffHcurlInvJac);
228 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
233 "This operator can be used only with element which is triangle");
243 const unsigned int nb_base_functions = data.
getDiffN(base).size2() / 6;
244 if (nb_base_functions) {
245 const unsigned int nb_gauss_pts = data.
getDiffN(base).size1();
247 diffHcurlInvJac.resize(nb_gauss_pts, nb_base_functions * 9,
false);
250 double *t_inv_diff_n_ptr = &*diffHcurlInvJac.data().begin();
252 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[
HVEC0_1],
261 auto t_inv_jac = getFTensor2FromMat<3, 3>(*invJacPtr);
262 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
263 for (
unsigned int bb = 0; bb != nb_base_functions; bb++) {
264 t_inv_diff_n(
i,
j) = t_diff_n(
i,
K) * t_inv_jac(
K,
j);
270 data.
getDiffN(base).swap(diffHcurlInvJac);
282 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
287 "This operator can be used only with element which is face");
293 const size_t nb_base_functions = data.
getN(base).size2() / 3;
294 if (nb_base_functions) {
296 const size_t nb_gauss_pts = data.
getN(base).size1();
300 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
301 for (
size_t bb = 0; bb != nb_base_functions; ++bb) {
303 const double a = t_n(0);
307 for (
auto n : {0, 1}) {
308 const double b = t_diff_n(0,
n);
309 t_diff_n(0,
n) = -t_diff_n(1,
n);
324 2>::OpSetCovariantPiolaTransformOnFace2DImpl(boost::shared_ptr<MatrixDouble>
327 invJacPtr(inv_jac_ptr) {}
333 const auto type_dim = moab::CN::Dimension(type);
334 if (type_dim != 1 && type_dim != 2)
345 auto &baseN = data.
getN(base);
346 auto &diffBaseN = data.
getDiffN(base);
348 int nb_dofs = baseN.size2() / 3;
349 int nb_gauss_pts = baseN.size1();
351 piolaN.resize(baseN.size1(), baseN.size2());
353 if (nb_dofs > 0 && nb_gauss_pts > 0) {
356 auto t_transformed_h_curl =
357 getFTensor1FromPtr<3, 3>(&*piolaN.data().begin());
358 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
359 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
360 for (
int ll = 0; ll != nb_dofs; ll++) {
361 t_transformed_h_curl(
i) = t_inv_jac(
j,
i) * t_h_curl(
j);
363 ++t_transformed_h_curl;
369 diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
370 if (diffBaseN.data().size() > 0) {
372 auto t_transformed_diff_h_curl =
374 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
375 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
376 for (
int ll = 0; ll != nb_dofs; ll++) {
377 t_transformed_diff_h_curl(
i,
k) =
378 t_inv_jac(
j,
i) * t_diff_h_curl(
j,
k);
380 ++t_transformed_diff_h_curl;
384 diffBaseN.swap(diffPiolaN);
397 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
408 const size_t nb_base_functions = data.
getN(base).size2() / 3;
409 if (nb_base_functions) {
411 const size_t nb_gauss_pts = data.
getN(base).size1();
412 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
413 if (data.
getN(base).size2() > 0) {
415 double *t_transformed_n_ptr = &*piolaN.data().begin();
418 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
419 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
420 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
423 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
424 t_transformed_n(
i) = t_jac(
i,
k) * t_n(
k) / det;
429 data.
getN(base).swap(piolaN);
432 piolaDiffN.resize(nb_gauss_pts, data.
getDiffN(base).size2(),
false);
433 if (data.
getDiffN(base).data().size() > 0) {
435 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
437 t_transformed_diff_n(t_transformed_diff_n_ptr,
438 &t_transformed_diff_n_ptr[
HVEC0_1],
440 &t_transformed_diff_n_ptr[
HVEC1_0],
441 &t_transformed_diff_n_ptr[
HVEC1_1],
443 &t_transformed_diff_n_ptr[
HVEC2_0],
444 &t_transformed_diff_n_ptr[
HVEC2_1]);
445 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
446 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
449 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
450 t_transformed_diff_n(
i,
k) = t_jac(
i,
j) * t_diff_n(
j,
k) / det;
452 ++t_transformed_diff_n;
455 data.
getDiffN(base).swap(piolaDiffN);
468 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
479 const size_t nb_base_functions = data.
getN(base).size2() / 3;
480 if (nb_base_functions) {
482 const size_t nb_gauss_pts = data.
getN(base).size1();
483 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
484 if (data.
getN(base).size2() > 0) {
486 double *t_transformed_n_ptr = &*piolaN.data().begin();
489 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
490 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
491 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
494 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
495 t_transformed_n(
i) = t_jac(
i,
j) * t_n(
j) / det;
500 data.
getN(base).swap(piolaN);
503 piolaDiffN.resize(nb_gauss_pts, data.
getDiffN(base).size2(),
false);
504 if (data.
getDiffN(base).size2() > 0) {
506 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
508 t_transformed_diff_n(t_transformed_diff_n_ptr,
509 &t_transformed_diff_n_ptr[
HVEC0_1],
511 &t_transformed_diff_n_ptr[
HVEC1_0],
512 &t_transformed_diff_n_ptr[
HVEC1_1],
514 &t_transformed_diff_n_ptr[
HVEC2_0],
515 &t_transformed_diff_n_ptr[
HVEC2_1]);
517 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
518 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
521 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
522 t_transformed_diff_n(
i,
K) = t_jac(
i,
j) * t_diff_n(
j,
K) / det;
524 ++t_transformed_diff_n;
527 data.
getDiffN(base).swap(piolaDiffN);
549 size_t nb_gauss_pts = data.
getN(base).size1();
550 size_t nb_dofs = data.
getN(base).size2() / 3;
557 if (edge_direction_at_gauss_pts.size1() == nb_gauss_pts) {
559 auto t_dir = getFTensor1TangentAtGaussPts<3>();
561 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
566 const auto l2 = t_normal(
i) * t_normal(
i);
567 for (
int ll = 0; ll != nb_dofs; ++ll) {
568 const double val = t_h_div(0);
569 const double a = val / l2;
570 t_h_div(
i) = t_normal(
i) *
a;
578 if (cc != nb_gauss_pts * nb_dofs)
590 if (type == MBVERTEX) {
593 double *coords_ptr = &*coords.data().begin();
595 const int nb_gauss_pts = data.
getN(
NOBASE).size1();
603 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
608 for (
int bb = 0; bb != 6; bb++) {
609 t_jac(
i,
j) += t_coords(
i) * t_diff_n(
j);
623 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
640 if (type == MBVERTEX) {
643 double *coords_ptr = &*coords.data().begin();
645 const int nb_gauss_pts = data.
getN(
NOBASE).size1();
647 invJac.resize(9, nb_gauss_pts,
false);
649 auto t_inv_jac = getFTensor2FromMat<3, 3>(
invJac);
655 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
660 for (
int bb = 0; bb != 6; bb++) {
661 t_jac(
i,
j) += t_coords(
i) * t_diff_n(
j);
687 if (data.
getN(base).size2() == 0)
690 const int nb_gauss_pts = data.
getN(base).size1();
700 auto t_inv_jac = getFTensor2FromMat<3, 3>(
invJac);
702 const int nb_dofs = data.
getN(base).size2();
703 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
704 for (
int bb = 0; bb != nb_dofs; bb++) {
705 t_inv_diff_n(
i) = t_diff_n(
j) * t_inv_jac(
j,
i);
724 if (type == MBVERTEX) {
727 double *coords_ptr = &*coords.data().begin();
730 double j00_f3, j01_f3, j10_f3, j11_f3;
731 for (
int gg = 0; gg < 1; gg++) {
734 j00_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[0], 2);
735 j01_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[1], 2);
736 j10_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[0], 2);
737 j11_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[1], 2);
739 double det_f3 = j00_f3 * j11_f3 - j01_f3 * j10_f3;
762 unsigned int nb_dofs = data.
getN(base).size2();
765 unsigned int nb_gauss_pts = data.
getN(base).size1();
766 diffNinvJac.resize(nb_gauss_pts, 2 * nb_dofs,
false);
769 if (type != MBVERTEX) {
770 if (nb_dofs != data.
getDiffN(base).size2() / 2) {
772 "data inconsistency nb_dofs != data.diffN.size2()/2 ( %u != "
774 nb_dofs, data.
getDiffN(base).size2());
783 for (
unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
784 for (
unsigned int dd = 0; dd < nb_dofs; dd++) {
785 cblas_dgemv(CblasRowMajor, CblasTrans, 2, 2, 1,
787 &data.
getDiffN(base)(gg, 2 * dd), 1, 0,
802 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
806 dataPtr(data_ptr), zeroType(
zero_type), zeroSide(zero_side) {
815 const auto nb_integration_points = getGaussPts().size2();
816 if (type == zeroType && side == zeroSide) {
817 dataPtr->resize(3, nb_integration_points,
false);
826 const auto nb_base_functions = data.
getN().size2() / 3;
828 auto t_data = getFTensor1FromMat<3>(*dataPtr);
829 for (
auto gg = 0; gg != nb_integration_points; ++gg) {
832 for (; bb != nb_dofs; ++bb) {
833 t_data(
k) += t_dof * (levi_civita(
j,
i,
k) * t_n_diff_hcurl(
i,
j));
837 for (; bb < nb_base_functions; ++bb)
846 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
850 dataPtr(data_ptr), zeroType(
zero_type), zeroSide(zero_side) {
859 const auto nb_integration_points = getGaussPts().size2();
860 if (type == zeroType && side == zeroSide) {
861 dataPtr->resize(2, nb_integration_points,
false);
871 const auto nb_base_functions = data.
getN().size2();
873 auto t_data = getFTensor1FromMat<2>(*dataPtr);
874 for (
auto gg = 0; gg != nb_integration_points; ++gg) {
877 for (; bb != nb_dofs; ++bb) {
878 t_data(
i) += t_dof * (levi_civita(
i,
j) * t_n_diff_hcurl(
j));
882 for (; bb < nb_base_functions; ++bb)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ USER_BASE
user implemented approximation base
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
FTensor::Tensor2< FTensor::PackPtr< double *, 4 >, 2, 2 > getFTensor2FromPtr< 2, 2 >(double *ptr)
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.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2FromPtr< 3, 2 >(double *ptr)
constexpr auto field_name
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
bool & doVertices
\deprectaed If false skip vertices
MatrixDouble & getTangentAtGaussPts()
get tangent vector to edge curve at integration points
static FTensor::Tensor1< double, 3 > tFaceOrientation
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()
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
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.
VectorDouble & getCoords()
get triangle coordinates
auto getFTensor0IntegrationWeight()
Get integration weights.
int getFEDim() const
Get dimension of finite element.
structure to get information form mofem into EntitiesFieldData
Calculate curl of vector field.
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.
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.
Transform local reference derivatives of shape function to global derivatives for face.
Transform local reference derivatives of shape functions to global derivatives.
OpSetInvJacToScalarBasesBasic(FieldSpace space, boost::shared_ptr< MatrixDouble > inv_jac_ptr)
VectorDouble & getCoords()
nodal coordinates
double getVolume() const
element volume (linear geometry)