 |
| v0.14.0
|
Go to the documentation of this file.
12 meshPositionsFieldName(
"MESH_NODE_POSITIONS"), coords(24), jAc(3, 3),
14 opContravariantPiolaTransform(elementMeasure, jAc),
15 opCovariantPiolaTransform(
invJac), opSetInvJacHdivAndHcurl(
invJac),
16 tJac(&jAc(0, 0), &jAc(0, 1), &jAc(0, 2), &jAc(1, 0), &jAc(1, 1),
17 &jAc(1, 2), &jAc(2, 0), &jAc(2, 1), &jAc(2, 2)),
34 auto get_rule_by_type = [&]() {
37 return getRule(order_row + 1, order_col + 1, order_data + 1);
39 return getRule(order_row, order_col, order_data);
43 const int rule = get_rule_by_type();
45 auto calc_base_for_tet = [&]() {
47 const size_t nb_gauss_pts =
gaussPts.size2();
50 base.resize(nb_gauss_pts, 4,
false);
51 diff_base.resize(nb_gauss_pts, 12,
false);
54 double *diff_shape_ptr = &*diff_base.data().begin();
55 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
56 for (
int nn = 0; nn != 4; ++nn) {
57 for (
int dd = 0;
dd != 3; ++
dd, ++diff_shape_ptr) {
65 auto calc_base_for_hex = [&]() {
67 const size_t nb_gauss_pts =
gaussPts.size2();
70 base.resize(nb_gauss_pts, 8,
false);
71 diff_base.resize(nb_gauss_pts, 24,
false);
72 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
112 auto set_integration_pts_for_hex = [&]() {
119 auto set_integration_pts_for_tet = [&]() {
130 gaussPts.resize(4, nb_gauss_pts,
false);
131 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[1], 4,
133 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[2], 4,
135 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[3], 4,
144 base.resize(nb_gauss_pts, 4,
false);
145 diff_base.resize(nb_gauss_pts, 12,
false);
146 double *shape_ptr = &*base.data().begin();
147 cblas_dcopy(4 * nb_gauss_pts,
QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
149 double *diff_shape_ptr = &*diff_base.data().begin();
150 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
151 for (
int nn = 0; nn != 4; ++nn) {
152 for (
int dd = 0;
dd != 3; ++
dd, ++diff_shape_ptr) {
168 CHKERR set_integration_pts_for_tet();
171 CHKERR set_integration_pts_for_hex();
172 CHKERR calc_base_for_hex();
176 "Element type not implemented: %d",
type);
181 const size_t nb_gauss_pts =
gaussPts.size2();
185 CHKERR calc_base_for_tet();
188 CHKERR calc_base_for_hex();
192 "Element type not implemented: %d",
type);
208 auto get_tet_t_diff_n = [&]() {
214 auto get_hex_t_diff_n = [&]() {
221 auto get_t_diff_n = [&]() {
223 return get_tet_t_diff_n();
224 return get_hex_t_diff_n();
227 auto t_diff_n = get_t_diff_n();
236 tJac(
i,
j) += t_coords(
i) * t_diff_n(
j);
263 double *shape_functions_ptr = &*shape_functions.data().begin();
264 const size_t nb_base_functions = shape_functions.size2();
265 const size_t nb_gauss_pts =
gaussPts.size2();
272 shape_functions_ptr);
273 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
276 for (
int bb = 0; bb != nb_base_functions; ++bb) {
277 t_coords_at_gauss_ptr(
i) += t_coords(
i) * t_shape_functions;
281 ++t_coords_at_gauss_ptr;
294 auto dim_type = CN::Dimension(
type);
296 auto get_data_on_ents = [&](
auto lower_dim,
auto space) {
301 for (
auto dd = dim_type;
dd >= lower_dim; --
dd) {
302 int nb_ents = moab::CN::NumSubEntities(
type,
dd);
303 for (
int ii = 0; ii != nb_ents; ++ii) {
304 auto sub_ent_type = moab::CN::SubEntityType(
type,
dd, ii);
306 auto &data_on_ent = data->dataOnEntities[sub_ent_type];
309 data->spacesOnEntities[sub_ent_type].set(space);
334 if ((data.
getDiffN(base).size1() == 4) &&
335 (data.
getDiffN(base).size2() == 3)) {
336 const size_t nb_gauss_pts =
gaussPts.size2();
337 const size_t nb_base_functions = 4;
338 new_diff_n.resize(nb_gauss_pts, 3 * nb_base_functions,
false);
339 double *new_diff_n_ptr = &*new_diff_n.data().begin();
341 new_diff_n_ptr, &new_diff_n_ptr[1], &new_diff_n_ptr[2]);
342 double *t_diff_n_ptr = &*data.
getDiffN(base).data().begin();
343 for (
unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
345 t_diff_n_ptr, &t_diff_n_ptr[1], &t_diff_n_ptr[2]);
346 for (
unsigned int bb = 0; bb != nb_base_functions; bb++) {
347 t_new_diff_n(
i) = t_diff_n(
i);
352 data.
getDiffN(base).resize(new_diff_n.size1(), new_diff_n.size2(),
354 data.
getDiffN(base).swap(new_diff_n);
362 std::array<std::bitset<LASTSPACE>, 4> spaces_by_dim{
364 std::bitset<LASTSPACE>(0),
365 std::bitset<LASTSPACE>(0),
366 std::bitset<LASTSPACE>(0),
367 std::bitset<LASTSPACE>(0)
376 spaces_by_dim[1].test(
HCURL) ||
377 spaces_by_dim[2].test(
HCURL) ||
378 spaces_by_dim[3].test(
HCURL)
385 if (spaces_by_dim[2].test(
HDIV) || spaces_by_dim[3].test(
HDIV)) {
391 for (
auto t : {MBTRI, MBTET}) {
394 d.getDiffN(base) /= 6;
402 if (spaces_by_dim[3].test(
L2)) {
414 "User operator and finite element do not work together");
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Data on single entity (This is passed as argument to DataOperator::doWork)
virtual MoFEMErrorCode setIntegrationPts()
Set integration points.
Calculate base functions on tetrahedral.
MatrixInt facesNodesOrder
order of face nodes on element
MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement()
Calculate Bernstein-Bezier base.
virtual MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
int getMaxRowOrder() const
Get max order of approximation for field in rows.
MoFEMErrorCode getFaceNodes(EntitiesFieldData &data) const
Get nodes on faces.
EntityType lastEvaluatedElementEntityType
Last evaluated type of element entity.
virtual MPI_Comm & get_comm() const =0
@ L2
field with C-1 continuity
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
OpSetCovariantPiolaTransform opCovariantPiolaTransform
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
int getMaxDataOrder() const
Get max order of approximation for data fields.
EntitiesFieldData & dataHdiv
OpSetInvJacHdivAndHcurl opSetInvJacHdivAndHcurl
virtual MoFEMErrorCode calculateCoordinatesAtGaussPts()
Calculate coordinate at integration points.
MoFEMErrorCode getEntityDataOrder(const EntityType type, const FieldSpace space, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
Get the entity data order.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
int getMaxColOrder() const
Get max order of approximation for field in columns.
MoFEMErrorCode calHierarchicalBaseFunctionsOnElement()
Calculate base functions.
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
MatrixInt facesNodes
nodes on finite element faces
double zeta
Viscous hardening.
Deprecated interface functions.
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.
EntitiesFieldData & dataL2
#define CHKERR
Inline error check.
virtual MoFEMErrorCode calculateVolumeAndJacobian()
Calculate element volume and Jacobian.
virtual moab::Interface & get_moab()=0
implementation of Data Operators for Forces and Sources
MatrixDouble coordsAtGaussPts
coordinated at gauss points
virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement()
Determine approximation space and order of base functions.
VolumeElementForcesAndSourcesCore(Interface &m_field)
MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr)
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
const EntityHandle * conn
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
OpSetInvJacH1 opSetInvJacH1
FTensor::Tensor2< double *, 3, 3 > tJac
constexpr double t
plate stiffness
Volume finite element base.
FTensor::Index< 'i', SPACE_DIM > i
static QUAD *const QUAD_3D_TABLE[]
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
structure to get information form mofem into EntitiesFieldData
#define QUAD_3D_TABLE_SIZE
EntitiesFieldData & dataH1
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
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)
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'j', 3 > j
EntitiesFieldData & dataHcurl
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ HCURL
field with continuous tangents
@ MOFEM_DATA_INCONSISTENCY
MatrixDouble gaussPts
Matrix of integration points.
OpSetContravariantPiolaTransform opContravariantPiolaTransform
FieldApproximationBase
approximation base
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
MoFEMErrorCode operator()()
function is run for every finite element
FTensor::Tensor2< double *, 3, 3 > tInvJac
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ HDIV
field with continuous normal traction
ForcesAndSourcesCore * ptrFE
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
virtual MoFEMErrorCode transformBaseFunctions()
Transform base functions based on geometric element Jacobian.
Calculate base functions on tetrahedral.
MoFEMErrorCode getEntitySense(const EntityType type, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
get sense (orientation) of entity