26 meshPositionsFieldName(
"MESH_NODE_POSITIONS") {}
43 auto get_ftensor_n_diff = [&]() {
51 &
m(0, 0), &
m(0, 1), &
m(0, 2));
56 const size_t nb_gauss_pts =
gaussPts.size2();
72 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
81 }
else if (
type == MBQUAD) {
89 const size_t nb_gauss_pts =
gaussPts.size2();
104 auto t_diff = get_ftensor_n_diff();
105 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
106 auto t_coords = get_ftensor_from_vec_3d(
coords);
107 for (
int nn = 0; nn !=
num_nodes; ++nn) {
108 t_t1(
i) += t_coords(
i) * t_diff(
N0);
109 t_t2(
i) += t_coords(
i) * t_diff(
N1);
121 "Element type not implemented");
138 auto calc_normal = [&](
const double *diff_ptr) {
149 diff_ptr, &diff_ptr[1]);
160 for (
int nn = 0; nn !=
num_nodes; ++nn) {
161 t_t1(
i) += t_coords(
i) * t_diff(
N0);
162 t_t2(
i) += t_coords(
i) * t_diff(
N1);
167 aRea = sqrt(t_normal(
i) * t_normal(
i));
171 const double *diff_ptr;
175 CHKERR calc_normal(diff_ptr);
181 CHKERR calc_normal(diff_ptr);
185 "Element type not implemented");
197 int rule =
getRule(order_row, order_col, order_data);
199 auto set_integration_pts_for_tri = [&]() {
210 gaussPts.resize(3, nb_gauss_pts,
false);
211 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[1], 3,
213 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[2], 3,
221 cblas_dcopy(3 * nb_gauss_pts,
QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
234 auto calc_base_for_quad = [&]() {
236 const size_t nb_gauss_pts =
gaussPts.size2();
239 base.resize(nb_gauss_pts, 4,
false);
240 diff_base.resize(nb_gauss_pts, 8,
false);
241 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
243 const double zeta =
gaussPts(1, gg);
265 CHKERR set_integration_pts_for_tri();
270 CHKERR calc_base_for_quad();
274 "Element type not implemented: %d",
type);
280 const size_t nb_gauss_pts =
gaussPts.size2();
286 base.resize(nb_gauss_pts, 3,
false);
287 diff_base.resize(3, 2,
false);
295 CHKERR calc_base_for_quad();
299 "Element type not implemented: %d",
type);
375 const size_t nb_nodes =
377 double *shape_functions =
379 const size_t nb_gauss_pts =
gaussPts.size2();
381 for (
int gg = 0; gg != nb_gauss_pts; ++gg)
382 for (
int dd = 0;
dd != 3; ++
dd)
384 nb_nodes, &shape_functions[nb_nodes * gg], 1, &
coords[
dd], 3);
394 "User operator and finite element do not work together");
441 return loopSide(fe_name, &fe_method, 3);
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#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 ...
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
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
#define QUAD_2D_TABLE_SIZE
static QUAD *const QUAD_2D_TABLE[]
FTensor::Index< 'm', 3 > m
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr)
MoFEMErrorCode loopSideVolumes(const string fe_name, VolumeElementForcesAndSourcesCoreOnSide &fe_method)
MatrixDouble tangentOneAtGaussPts
virtual MoFEMErrorCode calculateAreaAndNormal()
Calculate element area and normal of the face.
MoFEMErrorCode operator()()
function is run for every finite element
FaceElementForcesAndSourcesCore(Interface &m_field)
MatrixDouble normalsAtGaussPts
virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement()
Determine approximation space and order of base functions.
virtual MoFEMErrorCode calculateAreaAndNormalAtIntegrationPts()
Calculate element area and normal of the face at integration points.
virtual MoFEMErrorCode setIntegrationPts()
Set integration points.
const EntityHandle * conn
virtual MoFEMErrorCode calculateCoordinatesAtGaussPts()
Calculate coordinate at integration points.
MatrixDouble tangentTwoAtGaussPts
ForcesAndSourcesCore * ptrFE
structure to get information form mofem into EntitiesFieldData
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
int getMaxRowOrder() const
Get max order of approximation for field in rows.
EntitiesFieldData & dataHdiv
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
EntitiesFieldData & dataHcurl
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
EntitiesFieldData & dataH1
MatrixDouble coordsAtGaussPts
coordinated at gauss points
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
MoFEMErrorCode calHierarchicalBaseFunctionsOnElement()
Calculate base functions.
MatrixDouble gaussPts
Matrix of integration points.
MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement()
Calculate Bernstein-Bezier base.
int getMaxColOrder() const
Get max order of approximation for field in columns.
EntityType lastEvaluatedElementEntityType
Last evaluated type of element entity.
int getMaxDataOrder() const
Get max order of approximation for data fields.
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
EntitiesFieldData & dataL2
Calculate base functions on triangle.
Calculate base functions on triangle.
Base volume element used to integrate on skeleton.