11 FieldSpace space, boost::shared_ptr<MatrixDouble> inv_jac_ptr)
27 "This operator can be used only with element which faces");
30 return applyTransform<2, 2, 2, 2>(diff_n);
33 if (!(type == MBVERTEX && sPace ==
L2)) {
43 CHKERR apply_transform(*(
m.second));
48 CHKERR apply_transform(*ptr);
62 "This operator can be used only with element which face");
65 return applyTransform<2, 3, 3, 3>(diff_n);
68 if (!(type == MBVERTEX && sPace ==
L2)) {
77 CHKERR apply_transform(*(
m.second));
82 CHKERR apply_transform(*ptr);
96 "This operator can be used only with element which faces");
100 size_t nb_functions = diff2_n.size2() / 4;
102 size_t nb_gauss_pts = diff2_n.size1();
105 if (nb_gauss_pts != getGaussPts().size2())
107 "Wrong number of Gauss Pts");
110 diffNinvJac.resize(diff2_n.size1(), diff2_n.size2(),
false);
119 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
120 for (
size_t dd = 0; dd != nb_functions; ++dd) {
122 t_inv_jac(
k,
i) * (t_inv_jac(
l,
j) * t_diff2_n_ref(
k,
l));
128 diff2_n.swap(diffNinvJac);
133 if (!(type == MBVERTEX && sPace ==
L2)) {
138 data.
getN(base, BaseDerivatives::SecondDerivative));
141 auto error = [&](
auto &
m) {
145 "Operator do not work for Bernstein-Bezier. This "
155 CHKERR apply_transform(*(
m.second));
162 CHKERR apply_transform(*ptr);
175 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
180 "This operator can be used only with element which is triangle");
190 const unsigned int nb_base_functions = data.
getDiffN(base).size2() / 6;
191 if (nb_base_functions) {
192 const unsigned int nb_gauss_pts = data.
getDiffN(base).size1();
194 diffHcurlInvJac.resize(nb_gauss_pts, data.
getDiffN(base).size2(),
false);
197 double *t_inv_diff_n_ptr = &*diffHcurlInvJac.data().begin();
199 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[
HVEC0_1],
206 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
207 for (
unsigned int bb = 0; bb != nb_base_functions; bb++) {
208 t_inv_diff_n(
i,
j) = t_diff_n(
i,
k) * t_inv_jac(
k,
j);
214 data.
getDiffN(base).swap(diffHcurlInvJac);
226 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
231 "This operator can be used only with element which is triangle");
241 const unsigned int nb_base_functions = data.
getDiffN(base).size2() / 6;
242 if (nb_base_functions) {
243 const unsigned int nb_gauss_pts = data.
getDiffN(base).size1();
245 diffHcurlInvJac.resize(nb_gauss_pts, nb_base_functions * 9,
false);
248 double *t_inv_diff_n_ptr = &*diffHcurlInvJac.data().begin();
250 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[
HVEC0_1],
260 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
261 for (
unsigned int bb = 0; bb != nb_base_functions; bb++) {
262 t_inv_diff_n(
i,
j) = t_diff_n(
i,
K) * t_inv_jac(
K,
j);
268 data.
getDiffN(base).swap(diffHcurlInvJac);
279 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
284 "This operator can be used only with element which is face");
290 const size_t nb_base_functions = data.
getN(base).size2() / 3;
291 if (nb_base_functions) {
293 const size_t nb_gauss_pts = data.
getN(base).size1();
297 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
298 for (
size_t bb = 0; bb != nb_base_functions; ++bb) {
300 const double a = t_n(0);
304 for (
auto n : {0, 1}) {
305 const double b = t_diff_n(0,
n);
306 t_diff_n(0,
n) = -t_diff_n(1,
n);
321 2>::OpSetCovariantPiolaTransformOnFace2DImpl(boost::shared_ptr<MatrixDouble>
324 invJacPtr(inv_jac_ptr) {}
330 const auto type_dim = moab::CN::Dimension(type);
331 if (type_dim != 1 && type_dim != 2)
342 auto &baseN = data.
getN(base);
343 auto &diffBaseN = data.
getDiffN(base);
345 int nb_dofs = baseN.size2() / 3;
346 int nb_gauss_pts = baseN.size1();
348 piolaN.resize(baseN.size1(), baseN.size2());
350 if (nb_dofs > 0 && nb_gauss_pts > 0) {
353 auto t_transformed_h_curl =
356 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
357 for (
int ll = 0; ll != nb_dofs; ll++) {
358 t_transformed_h_curl(
i) = t_inv_jac(
j,
i) * t_h_curl(
j);
360 ++t_transformed_h_curl;
366 diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
367 if (diffBaseN.data().size() > 0) {
369 auto t_transformed_diff_h_curl =
372 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
373 for (
int ll = 0; ll != nb_dofs; ll++) {
374 t_transformed_diff_h_curl(
i,
k) =
375 t_inv_jac(
j,
i) * t_diff_h_curl(
j,
k);
377 ++t_transformed_diff_h_curl;
381 diffBaseN.swap(diffPiolaN);
394 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
405 const size_t nb_base_functions = data.
getN(base).size2() / 3;
406 if (nb_base_functions) {
408 const size_t nb_gauss_pts = data.
getN(base).size1();
409 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
410 if (data.
getN(base).size2() > 0) {
412 double *t_transformed_n_ptr = &*piolaN.data().begin();
415 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
417 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
420 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
421 t_transformed_n(
i) = t_jac(
i,
k) * t_n(
k) / det;
426 data.
getN(base).swap(piolaN);
429 piolaDiffN.resize(nb_gauss_pts, data.
getDiffN(base).size2(),
false);
430 if (data.
getDiffN(base).data().size() > 0) {
432 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
434 t_transformed_diff_n(t_transformed_diff_n_ptr,
435 &t_transformed_diff_n_ptr[
HVEC0_1],
437 &t_transformed_diff_n_ptr[
HVEC1_0],
438 &t_transformed_diff_n_ptr[
HVEC1_1],
440 &t_transformed_diff_n_ptr[
HVEC2_0],
441 &t_transformed_diff_n_ptr[
HVEC2_1]);
443 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
446 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
447 t_transformed_diff_n(
i,
k) = t_jac(
i,
j) * t_diff_n(
j,
k) / det;
449 ++t_transformed_diff_n;
452 data.
getDiffN(base).swap(piolaDiffN);
465 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
476 const size_t nb_base_functions = data.
getN(base).size2() / 3;
477 if (nb_base_functions) {
479 const size_t nb_gauss_pts = data.
getN(base).size1();
480 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
481 if (data.
getN(base).size2() > 0) {
483 double *t_transformed_n_ptr = &*piolaN.data().begin();
486 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
488 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
491 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
492 t_transformed_n(
i) = t_jac(
i,
j) * t_n(
j) / det;
497 data.
getN(base).swap(piolaN);
500 piolaDiffN.resize(nb_gauss_pts, data.
getDiffN(base).size2(),
false);
501 if (data.
getDiffN(base).size2() > 0) {
503 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
505 t_transformed_diff_n(t_transformed_diff_n_ptr,
506 &t_transformed_diff_n_ptr[
HVEC0_1],
508 &t_transformed_diff_n_ptr[
HVEC1_0],
509 &t_transformed_diff_n_ptr[
HVEC1_1],
511 &t_transformed_diff_n_ptr[
HVEC2_0],
512 &t_transformed_diff_n_ptr[
HVEC2_1]);
515 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
518 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
519 t_transformed_diff_n(
i,
K) = t_jac(
i,
j) * t_diff_n(
j,
K) / det;
521 ++t_transformed_diff_n;
524 data.
getDiffN(base).swap(piolaDiffN);
546 size_t nb_gauss_pts = data.
getN(base).size1();
547 size_t nb_dofs = data.
getN(base).size2() / 3;
554 if (edge_direction_at_gauss_pts.size1() == nb_gauss_pts) {
558 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
563 const auto l2 = t_normal(
i) * t_normal(
i);
564 for (
int ll = 0; ll != nb_dofs; ++ll) {
565 const double val = t_h_div(0);
566 const double a = val / l2;
567 t_h_div(
i) = t_normal(
i) *
a;
575 if (cc != nb_gauss_pts * nb_dofs)
587 if (type == MBVERTEX) {
590 double *coords_ptr = &*coords.data().begin();
592 const int nb_gauss_pts = data.
getN(
NOBASE).size1();
600 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
605 for (
int bb = 0; bb != 6; bb++) {
606 t_jac(
i,
j) += t_coords(
i) * t_diff_n(
j);
620 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
637 if (type == MBVERTEX) {
640 double *coords_ptr = &*coords.data().begin();
642 const int nb_gauss_pts = data.
getN(
NOBASE).size1();
644 invJac.resize(9, nb_gauss_pts,
false);
652 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
657 for (
int bb = 0; bb != 6; bb++) {
658 t_jac(
i,
j) += t_coords(
i) * t_diff_n(
j);
684 if (data.
getN(base).size2() == 0)
687 const int nb_gauss_pts = data.
getN(base).size1();
699 const int nb_dofs = data.
getN(base).size2();
700 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
701 for (
int bb = 0; bb != nb_dofs; bb++) {
702 t_inv_diff_n(
i) = t_diff_n(
j) * t_inv_jac(
j,
i);
721 if (type == MBVERTEX) {
724 double *coords_ptr = &*coords.data().begin();
727 double j00_f3, j01_f3, j10_f3, j11_f3;
728 for (
int gg = 0; gg < 1; gg++) {
731 j00_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[0], 2);
732 j01_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[1], 2);
733 j10_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[0], 2);
734 j11_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[1], 2);
736 double det_f3 = j00_f3 * j11_f3 - j01_f3 * j10_f3;
759 auto nb_dofs = data.
getN(base).size2();
762 unsigned int nb_gauss_pts = data.
getN(base).size1();
763 diffNinvJac.resize(nb_gauss_pts, 2 * nb_dofs,
false);
766 if (type != MBVERTEX) {
767 if (nb_dofs != data.
getDiffN(base).size2() / 2) {
769 "data inconsistency nb_dofs != data.diffN.size2()/2 ( %ld != "
771 nb_dofs, data.
getDiffN(base).size2());
780 for (
unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
781 for (
unsigned int dd = 0; dd < nb_dofs; dd++) {
782 cblas_dgemv(CblasRowMajor, CblasTrans, 2, 2, 1,
784 &data.
getDiffN(base)(gg, 2 * dd), 1, 0,
799 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
800 const EntityType zero_type,
const int zero_side)
803 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
812 const auto nb_integration_points = getGaussPts().size2();
813 if (type == zeroType && side == zeroSide) {
814 dataPtr->resize(3, nb_integration_points,
false);
823 const auto nb_base_functions = data.
getN().size2() / 3;
826 for (
auto gg = 0; gg != nb_integration_points; ++gg) {
829 for (; bb != nb_dofs; ++bb) {
830 t_data(
k) += t_dof * (levi_civita(
j,
i,
k) * t_n_diff_hcurl(
i,
j));
834 for (; bb < nb_base_functions; ++bb)
843 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
844 const EntityType zero_type,
const int zero_side)
847 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
856 const auto nb_integration_points = getGaussPts().size2();
857 if (type == zeroType && side == zeroSide) {
858 dataPtr->resize(2, nb_integration_points,
false);
868 const auto nb_base_functions = data.
getN().size2();
871 for (
auto gg = 0; gg != nb_integration_points; ++gg) {
874 for (; bb != nb_dofs; ++bb) {
875 t_data(
i) += t_dof * (levi_civita(
i,
j) * t_n_diff_hcurl(
j));
879 for (; bb < nb_base_functions; ++bb)
888 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
889 const EntityType zero_type,
const int zero_side)
892 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
901 const auto nb_integration_points = getGaussPts().size2();
902 if (type == zeroType && side == zeroSide) {
903 dataPtr->resize(3, nb_integration_points,
false);
914 const auto nb_base_functions = data.
getN().size2();
917 auto t_normal = getFTensor1NormalsAtGaussPts();
918 for (
auto gg = 0; gg != nb_integration_points; ++gg) {
920 auto nrm = t_normal.l2();
922 for (; bb != nb_dofs; ++bb) {
923 t_data(
k) += (t_dof / nrm) *
924 ((levi_civita(
i,
j,
k) * t_normal(
i)) * t_n_diff_hcurl(
j));
928 for (; bb < nb_base_functions; ++bb)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
@ 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< 'i', SPACE_DIM > i
const double n
refractive index of diffusive medium
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
auto getFTensor2FromPtr(double *ptr)
Make Tensor2 from pointer.
static auto determinantTensor2by2(T &t)
Calculate the determinant of a 2x2 matrix or a tensor of rank 2.
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2HVecFromPtr< 3, 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.
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
constexpr auto field_name
FTensor::Index< 'm', 3 > m
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
FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, DIM > getFTensor1TangentAtGaussPts()
Get tangent vector 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 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.
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)