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)
337 const auto nb_gauss_pts = getGaussPts().size2();
347 auto &baseN = data.
getN(base);
348 auto &diffBaseN = data.
getDiffN(base);
350 int nb_dofs = baseN.size2() / 3;
351 int nb_gauss_pts = baseN.size1();
353 piolaN.resize(baseN.size1(), baseN.size2());
355 if (nb_dofs > 0 && nb_gauss_pts > 0) {
358 auto t_transformed_h_curl =
359 getFTensor1FromPtr<3, 3>(&*piolaN.data().begin());
360 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
361 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
362 for (
int ll = 0; ll != nb_dofs; ll++) {
363 t_transformed_h_curl(
i) = t_inv_jac(
j,
i) * t_h_curl(
j);
365 ++t_transformed_h_curl;
371 diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
372 if (diffBaseN.data().size() > 0) {
374 auto t_transformed_diff_h_curl =
376 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
377 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
378 for (
int ll = 0; ll != nb_dofs; ll++) {
379 t_transformed_diff_h_curl(
i,
k) =
380 t_inv_jac(
j,
i) * t_diff_h_curl(
j,
k);
382 ++t_transformed_diff_h_curl;
386 diffBaseN.swap(diffPiolaN);
399 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
410 const size_t nb_base_functions = data.
getN(base).size2() / 3;
411 if (nb_base_functions) {
413 const size_t nb_gauss_pts = data.
getN(base).size1();
414 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
415 if (data.
getN(base).size2() > 0) {
417 double *t_transformed_n_ptr = &*piolaN.data().begin();
420 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
421 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
422 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
425 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
426 t_transformed_n(
i) = t_jac(
i,
k) * t_n(
k) / det;
431 data.
getN(base).swap(piolaN);
434 piolaDiffN.resize(nb_gauss_pts, data.
getDiffN(base).size2(),
false);
435 if (data.
getDiffN(base).data().size() > 0) {
437 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
439 t_transformed_diff_n(t_transformed_diff_n_ptr,
440 &t_transformed_diff_n_ptr[
HVEC0_1],
442 &t_transformed_diff_n_ptr[
HVEC1_0],
443 &t_transformed_diff_n_ptr[
HVEC1_1],
445 &t_transformed_diff_n_ptr[
HVEC2_0],
446 &t_transformed_diff_n_ptr[
HVEC2_1]);
447 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
448 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
451 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
452 t_transformed_diff_n(
i,
k) = t_jac(
i,
j) * t_diff_n(
j,
k) / det;
454 ++t_transformed_diff_n;
457 data.
getDiffN(base).swap(piolaDiffN);
470 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
481 const size_t nb_base_functions = data.
getN(base).size2() / 3;
482 if (nb_base_functions) {
484 const size_t nb_gauss_pts = data.
getN(base).size1();
485 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
486 if (data.
getN(base).size2() > 0) {
488 double *t_transformed_n_ptr = &*piolaN.data().begin();
491 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
492 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
493 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
496 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
497 t_transformed_n(
i) = t_jac(
i,
j) * t_n(
j) / det;
502 data.
getN(base).swap(piolaN);
505 piolaDiffN.resize(nb_gauss_pts, data.
getDiffN(base).size2(),
false);
506 if (data.
getDiffN(base).size2() > 0) {
508 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
510 t_transformed_diff_n(t_transformed_diff_n_ptr,
511 &t_transformed_diff_n_ptr[
HVEC0_1],
513 &t_transformed_diff_n_ptr[
HVEC1_0],
514 &t_transformed_diff_n_ptr[
HVEC1_1],
516 &t_transformed_diff_n_ptr[
HVEC2_0],
517 &t_transformed_diff_n_ptr[
HVEC2_1]);
519 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
520 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
523 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
524 t_transformed_diff_n(
i,
K) = t_jac(
i,
j) * t_diff_n(
j,
K) / det;
526 ++t_transformed_diff_n;
529 data.
getDiffN(base).swap(piolaDiffN);
551 size_t nb_gauss_pts = data.
getN(base).size1();
552 size_t nb_dofs = data.
getN(base).size2() / 3;
559 if (edge_direction_at_gauss_pts.size1() == nb_gauss_pts) {
561 auto t_dir = getFTensor1TangentAtGaussPts<3>();
563 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
568 const auto l2 = t_normal(
i) * t_normal(
i);
569 for (
int ll = 0; ll != nb_dofs; ++ll) {
570 const double val = t_h_div(0);
571 const double a = val / l2;
572 t_h_div(
i) = t_normal(
i) *
a;
580 if (cc != nb_gauss_pts * nb_dofs)
592 if (type == MBVERTEX) {
595 double *coords_ptr = &*coords.data().begin();
597 const int nb_gauss_pts = data.
getN(
NOBASE).size1();
605 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
610 for (
int bb = 0; bb != 6; bb++) {
611 t_jac(
i,
j) += t_coords(
i) * t_diff_n(
j);
625 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
642 if (type == MBVERTEX) {
645 double *coords_ptr = &*coords.data().begin();
647 const int nb_gauss_pts = data.
getN(
NOBASE).size1();
649 invJac.resize(9, nb_gauss_pts,
false);
651 auto t_inv_jac = getFTensor2FromMat<3, 3>(
invJac);
658 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
663 for (
int bb = 0; bb != 6; bb++) {
664 t_jac(
i,
j) += t_coords(
i) * t_diff_n(
j);
690 if (data.
getN(base).size2() == 0)
693 const int nb_gauss_pts = data.
getN(base).size1();
703 auto t_inv_jac = getFTensor2FromMat<3, 3>(
invJac);
705 const int nb_dofs = data.
getN(base).size2();
706 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
707 for (
int bb = 0; bb != nb_dofs; bb++) {
708 t_inv_diff_n(
i) = t_diff_n(
j) * t_inv_jac(
j,
i);
727 if (type == MBVERTEX) {
730 double *coords_ptr = &*coords.data().begin();
733 double j00_f3, j01_f3, j10_f3, j11_f3;
734 for (
int gg = 0; gg < 1; gg++) {
737 j00_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[0], 2);
738 j01_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[1], 2);
739 j10_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[0], 2);
740 j11_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[1], 2);
742 double det_f3 = j00_f3 * j11_f3 - j01_f3 * j10_f3;
765 unsigned int nb_dofs = data.
getN(base).size2();
768 unsigned int nb_gauss_pts = data.
getN(base).size1();
769 diffNinvJac.resize(nb_gauss_pts, 2 * nb_dofs,
false);
772 if (type != MBVERTEX) {
773 if (nb_dofs != data.
getDiffN(base).size2() / 2) {
775 "data inconsistency nb_dofs != data.diffN.size2()/2 ( %u != "
777 nb_dofs, data.
getDiffN(base).size2());
786 for (
unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
787 for (
unsigned int dd = 0; dd < nb_dofs; dd++) {
788 cblas_dgemv(CblasRowMajor, CblasTrans, 2, 2, 1,
790 &data.
getDiffN(base)(gg, 2 * dd), 1, 0,
805 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
809 dataPtr(data_ptr), zeroType(
zero_type), zeroSide(zero_side) {
818 const auto nb_integration_points = getGaussPts().size2();
819 if (type == zeroType && side == zeroSide) {
820 dataPtr->resize(3, nb_integration_points,
false);
829 const auto nb_base_functions = data.
getN().size2() / 3;
831 auto t_data = getFTensor1FromMat<3>(*dataPtr);
832 for (
auto gg = 0; gg != nb_integration_points; ++gg) {
835 for (; bb != nb_dofs; ++bb) {
836 t_data(
k) += t_dof * (levi_civita(
j,
i,
k) * t_n_diff_hcurl(
i,
j));
840 for (; bb < nb_base_functions; ++bb)
849 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
853 dataPtr(data_ptr), zeroType(
zero_type), zeroSide(zero_side) {
862 const auto nb_integration_points = getGaussPts().size2();
863 if (type == zeroType && side == zeroSide) {
864 dataPtr->resize(2, nb_integration_points,
false);
874 const auto nb_base_functions = data.
getN().size2();
876 auto t_data = getFTensor1FromMat<2>(*dataPtr);
877 for (
auto gg = 0; gg != nb_integration_points; ++gg) {
880 for (; bb != nb_dofs; ++bb) {
881 t_data(
i) += t_dof * (levi_civita(
i,
j) * t_n_diff_hcurl(
j));
885 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)