24 boost::shared_ptr<MatrixDouble> jac_ptr)
28 for (
auto t = MBEDGE;
t != MBMAXTYPE; ++
t)
40 const auto nb_base_functions = data.
getN(
NOBASE).size2();
41 if (nb_base_functions) {
48 if (nb_gauss_pts != data.
getDiffN().size1())
50 "Inconsistent number base functions and gauss points");
51 if (nb_base_functions != data.
getDiffN().size2() / 3)
53 "Inconsistent number of base functions");
54 if (coords.size() != 3 * nb_base_functions)
57 "Number of vertex coordinates and base functions is inconsistent "
59 coords.size() / 3, nb_base_functions);
62 jacPtr->resize(9, nb_gauss_pts,
false);
64 auto t_jac = getFTensor2FromMat<3, 3>(*
jacPtr);
71 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
73 &coords[0], &coords[1], &coords[2]);
74 for (
size_t bb = 0; bb != nb_base_functions; ++bb) {
75 t_jac(
i,
j) += t_coords(
i) * (t_vol_inv_jac(
k,
j) * t_diff_base(
k));
97 const auto nb_integration_pts = data.
getN().size1();
98 const auto nb_base_functions = data.
getN().size2();
99 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
102 for (; bb != nb_dofs; ++bb) {
103 t_coords(
i) += t_base * t_dof(
i);
107 for (; bb < nb_base_functions; ++bb)
124 return applyTransform<3, 3, 3, 3>(diff_n);
136 CHKERR transform_base(*(
m.second));
141 CHKERR transform_base(*ptr);
165 data.
getDiffN(base).size2(),
false);
167 unsigned int nb_gauss_pts = data.
getDiffN(base).size1();
168 unsigned int nb_base_functions = data.
getDiffN(base).size2() / 9;
169 if (nb_base_functions == 0)
175 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[
HVEC0_1],
180 auto t_inv_jac = getFTensor2FromMat<3, 3>(*
invJacPtr);
181 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
182 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
183 t_inv_diff_n(
i,
j) = t_inv_jac(
k,
j) * t_diff_n(
i,
k);
209 for (
size_t gg = 0; gg != nb_int_pts; ++gg) {
210 t_w *= sqrt(t_normal(
i) * t_normal(
i)) /
a;
216 "Number of rows in getNormalsAtGaussPts should be equal to "
217 "number of integration points, but is not, i.e. %d != %d",
228 const auto nb_integration_pts =
detPtr->size();
233 "Inconsistent number of data points");
238 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
258 "Pointer for detPtr not allocated");
261 "Pointer for jacPtr not allocated");
268 auto nb_gauss_pts = data.
getN(base).size1();
269 auto nb_base_functions = data.
getN(base).size2() / 3;
272 if (data.
getDiffN(base).size1() != nb_gauss_pts)
274 "Wrong number integration points");
276 if (data.
getDiffN(base).size2() / 9 != nb_base_functions)
278 "Wrong number base functions %d != %d",
279 data.
getDiffN(base).size2(), nb_base_functions);
282 if (nb_gauss_pts && nb_base_functions) {
284 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
288 double *t_transformed_n_ptr = &*
piolaN.data().begin();
291 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
293 double *t_transformed_diff_n_ptr = &*
piolaDiffN.data().begin();
295 t_transformed_diff_n(t_transformed_diff_n_ptr,
296 &t_transformed_diff_n_ptr[
HVEC0_1],
297 &t_transformed_diff_n_ptr[
HVEC0_2],
298 &t_transformed_diff_n_ptr[
HVEC1_0],
299 &t_transformed_diff_n_ptr[
HVEC1_1],
300 &t_transformed_diff_n_ptr[
HVEC1_2],
301 &t_transformed_diff_n_ptr[
HVEC2_0],
302 &t_transformed_diff_n_ptr[
HVEC2_1],
303 &t_transformed_diff_n_ptr[
HVEC2_2]);
306 auto t_jac = getFTensor2FromMat<3, 3>(*
jacPtr);
308 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
309 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
310 const double a = 1. / t_det;
311 t_transformed_n(
i) =
a * t_jac(
i,
k) * t_n(
k);
312 t_transformed_diff_n(
i,
k) =
a * t_jac(
i,
j) * t_diff_n(
j,
k);
316 ++t_transformed_diff_n;
342 unsigned int nb_gauss_pts = data.
getN(base).size1();
343 unsigned int nb_base_functions = data.
getN(base).size2() / 3;
344 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
348 double *t_transformed_n_ptr = &*
piolaN.data().begin();
351 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
353 double *t_transformed_diff_n_ptr = &*
piolaDiffN.data().begin();
355 t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[
HVEC0_1],
356 &t_transformed_diff_n_ptr[
HVEC0_2], &t_transformed_diff_n_ptr[
HVEC1_0],
357 &t_transformed_diff_n_ptr[
HVEC1_1], &t_transformed_diff_n_ptr[
HVEC1_2],
358 &t_transformed_diff_n_ptr[
HVEC2_0], &t_transformed_diff_n_ptr[
HVEC2_1],
359 &t_transformed_diff_n_ptr[
HVEC2_2]);
361 auto t_inv_jac = getFTensor2FromMat<3, 3>(*
jacInvPtr);
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) = t_inv_jac(
k,
i) * t_n(
k);
366 t_transformed_diff_n(
i,
k) = t_inv_jac(
j,
i) * t_diff_n(
j,
k);
370 ++t_transformed_diff_n;
383 boost::shared_ptr<MatrixDouble> jac_ptr)
392 const auto nb_gauss_pts = getGaussPts().size2();
394 jacPtr->resize(4, nb_gauss_pts,
false);
395 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
400 auto t_t1 = getFTensor1Tangent1AtGaussPts();
401 auto t_t2 = getFTensor1Tangent2AtGaussPts();
405 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
406 t_jac(
i,
N0) = t_t1(
i);
407 t_jac(
i,
N1) = t_t2(
i);
424 size_t nb_gauss_pts = getGaussPts().size2();
425 jacPtr->resize(9, nb_gauss_pts,
false);
427 auto t_t1 = getFTensor1Tangent1AtGaussPts();
428 auto t_t2 = getFTensor1Tangent2AtGaussPts();
429 auto t_normal = getFTensor1NormalsAtGaussPts();
435 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
436 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
438 t_jac(
i,
N0) = t_t1(
i);
439 t_jac(
i,
N1) = t_t2(
i);
441 const double l = sqrt(t_normal(
j) * t_normal(
j));
443 t_jac(
i,
N2) = t_normal(
i) /
l;
471 &
m(0, 0), &
m(0, 1), &
m(0, 2));
477 int nb_gauss_pts = data.
getN().size1();
480 tangent1_at_gauss_pts.resize(nb_gauss_pts, 3,
false);
481 tangent2_at_gauss_pts.resize(nb_gauss_pts, 3,
false);
485 tangent1_at_gauss_pts.clear();
486 tangent2_at_gauss_pts.clear();
493 if (2 * data.
getN().size2() != data.
getDiffN().size2()) {
495 "data inconsistency");
497 if (nb_dofs % 3 != 0) {
499 "data inconsistency");
503 if (nb_dofs > 3 * data.
getN().size2()) {
505 for (; nn != nb_dofs; nn++) {
509 if (nn > 3 * data.
getN().size2()) {
511 "data inconsistency for base %s",
519 const int nb_base_functions = data.
getN().size2();
522 auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
523 auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
524 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
527 for (; bb != nb_dofs / 3; ++bb) {
528 t_t1(
i) += t_data(
i) * t_diff_base(
N0);
529 t_t2(
i) += t_data(
i) * t_diff_base(
N1);
534 for (; bb != nb_base_functions; ++bb) {
547 if (moab::CN::Dimension(
type) == 2) {
553 const auto nb_gauss_pts = tangent1_at_gauss_pts.size1();
554 normal_at_gauss_pts.resize(nb_gauss_pts, 3,
false);
556 auto t_normal = get_ftensor1(normal_at_gauss_pts);
557 auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
558 auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
559 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
575 if (moab::CN::Dimension(
type) != 2)
578 auto get_normals_ptr = [&]() {
585 auto apply_transform_nonlinear_geometry = [&](
auto base,
auto nb_gauss_pts,
586 auto nb_base_functions) {
589 auto ptr = get_normals_ptr();
591 &ptr[0], &ptr[1], &ptr[2]);
594 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
595 const auto l2 = t_normal(
i) * t_normal(
i);
596 for (
int bb = 0; bb != nb_base_functions; ++bb) {
597 const auto v = t_base(0);
598 t_base(
i) = (
v / l2) * t_normal(
i);
609 const auto &base_functions = data.
getN(base);
610 const auto nb_gauss_pts = base_functions.size1();
614 const auto nb_base_functions = base_functions.size2() / 3;
615 CHKERR apply_transform_nonlinear_geometry(base, nb_gauss_pts,
627 const auto type_dim = moab::CN::Dimension(
type);
628 if (type_dim != 1 && type_dim != 2)
641 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
643 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
645 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
651 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
653 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
655 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
663 auto &baseN = data.
getN(base);
664 auto &diffBaseN = data.
getDiffN(base);
666 int nb_dofs = baseN.size2() / 3;
667 int nb_gauss_pts = baseN.size1();
669 piolaN.resize(baseN.size1(), baseN.size2());
670 diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
672 if (nb_dofs > 0 && nb_gauss_pts > 0) {
681 t_transformed_diff_h_curl(
692 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
695 for (
int ll = 0; ll != nb_dofs; ll++) {
696 t_transformed_h_curl(
i) = t_inv_m(
j,
i) * t_h_curl(
j);
697 t_transformed_diff_h_curl(
i,
k) = t_inv_m(
j,
i) * t_diff_h_curl(
j,
k);
699 ++t_transformed_h_curl;
701 ++t_transformed_diff_h_curl;
708 if (cc != nb_gauss_pts * nb_dofs)
733 "Wrong number of integration points %d!=%d",
736 auto get_base_at_pts = [&](
auto base) {
743 auto get_tangent_ptr = [&]() {
751 auto get_tangent_at_pts = [&]() {
752 auto ptr = get_tangent_ptr();
754 &ptr[0], &ptr[1], &ptr[2]);
758 auto calculate_squared_edge_length = [&]() {
759 std::vector<double> l1;
761 l1.resize(nb_gauss_pts);
762 auto t_m_at_pts = get_tangent_at_pts();
763 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
764 l1[gg] = t_m_at_pts(
i) * t_m_at_pts(
i);
771 auto l1 = calculate_squared_edge_length();
776 const auto nb_dofs = data.
getN(base).size2() / 3;
778 if (nb_gauss_pts && nb_dofs) {
779 auto t_h_curl = get_base_at_pts(base);
781 auto t_m_at_pts = get_tangent_at_pts();
782 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
783 const double l0 = l1[gg];
784 for (
int ll = 0; ll != nb_dofs; ++ll) {
785 const double val = t_h_curl(0);
786 const double a = val / l0;
787 t_h_curl(
i) = t_m_at_pts(
i) *
a;
795 if (cc != nb_gauss_pts * nb_dofs)
820 auto &tangent = get_tangent();
822 int nb_gauss_pts = data.
getN().size1();
823 tangent.resize(nb_gauss_pts, 3,
false);
825 int nb_approx_fun = data.
getN().size2();
826 double *diff = &*data.
getDiffN().data().begin();
830 tangent.resize(nb_gauss_pts, 3,
false);
834 for (
int dd = 0;
dd != 3;
dd++) {
835 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
836 tangent(gg,
dd) = cblas_ddot(2, diff, 1, dofs[
dd], 3);
844 "Approximated field should be rank 3, i.e. vector in 3d space");
847 for (
int dd = 0;
dd != 3;
dd++) {
848 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
850 cblas_ddot(nb_dofs / 3, &diff[gg * nb_approx_fun], 1, dofs[
dd], 3);
856 "This operator can calculate tangent vector only on edge");
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static MoFEMErrorCode get_jac(EntitiesFieldData::EntData &col_data, int gg, MatrixDouble &jac_stress, MatrixDouble &jac)
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#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 ...
@ MOFEM_DATA_INCONSISTENCY
static const char *const ApproximationBaseNames[]
#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
double v
phase velocity of light in medium (cm/ns)
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
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
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 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 double t
plate stiffness
constexpr auto field_name
FTensor::Index< 'm', 3 > m
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
MatrixDouble & getTangetAtGaussPts()
get tangent vector to edge curve at integration points
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::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
double getMeasure()
get measure of element
MatrixDouble & getTangent2AtGaussPts()
if higher order geometry return tangent vector to triangle at Gauss pts.
auto getFTensor1NormalsAtGaussPts()
get normal at integration points
MatrixDouble & getTangent1AtGaussPts()
if higher order geometry return tangent vector to triangle at Gauss pts.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
auto getFTensor0IntegrationWeight()
Get integration weights.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
int getFEDim() const
Get dimension of finite element.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate jacobian for face element.
boost::shared_ptr< MatrixDouble > jacPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHOJacForVolume(boost::shared_ptr< MatrixDouble > jac_ptr)
OpGetHONormalsOnFace(std::string field_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > tangentsAtPts
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.
boost::shared_ptr< MatrixDouble > invJacPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MatrixDouble diffHdivInvJac
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< VectorDouble > detPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
VectorDouble & getCoords()
nodal coordinates
FTensor::Tensor2< double *, 3, 3 > & getInvJac()
get element inverse Jacobian
Volume finite element base.