9#ifndef __CONVECTIVE_MASS_ELEMENT_HPP
10#define __CONVECTIVE_MASS_ELEMENT_HPP
13#error "MoFEM need to be compiled with ADOL-C"
94 auto add_ops = [&](
auto &fe) {
121 std::map<int, BlockData>
144 std::vector<VectorDouble>
valT;
146 std::vector<MatrixDouble>
jacT;
153 :
public VolumeElementForcesAndSourcesCore::UserDataOperator {
160 std::vector<VectorDouble> &values_at_gauss_pts,
161 std::vector<MatrixDouble> &gardient_at_gauss_pts);
167 EntitiesFieldData::EntData &data);
178 :
public VolumeElementForcesAndSourcesCore::UserDataOperator,
192 boost::ptr_vector<MethodForForceScaling> &methods_op,
193 int tag,
bool linear =
false);
204 EntitiesFieldData::EntData &row_data);
207 struct OpMassRhs :
public VolumeElementForcesAndSourcesCore::UserDataOperator,
219 EntitiesFieldData::EntData &row_data);
223 :
public VolumeElementForcesAndSourcesCore::UserDataOperator,
232 Range *forcesonlyonentities_ptr = NULL);
241 EntitiesFieldData::EntData &row_data,
242 EntitiesFieldData::EntData &col_data);
261 struct OpEnergy :
public VolumeElementForcesAndSourcesCore::UserDataOperator,
276 EntitiesFieldData::EntData &row_data);
280 :
public VolumeElementForcesAndSourcesCore::UserDataOperator,
289 CommonData &common_data,
int tag,
bool jacobian =
true);
298 EntitiesFieldData::EntData &row_data);
302 :
public VolumeElementForcesAndSourcesCore::UserDataOperator,
314 EntitiesFieldData::EntData &row_data);
348 :
public VolumeElementForcesAndSourcesCore::UserDataOperator,
360 bool jacobian =
true);
367 EntitiesFieldData::EntData &row_data);
371 :
public VolumeElementForcesAndSourcesCore::UserDataOperator,
381 Range *forcesonlyonentities_ptr);
386 EntitiesFieldData::EntData &row_data);
395 Range *forcesonlyonentities_ptr);
408 Range *forcesonlyonentities_ptr);
421 Range *forcesonlyonentities_ptr);
443 const std::string velocity_field,
444 const std::string spatial_position_field);
461 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr);
464 string element_name,
string velocity_field_name,
465 string spatial_position_field_name,
466 string material_position_field_name =
"MESH_NODE_POSITIONS",
470 string element_name,
string velocity_field_name,
471 string spatial_position_field_name,
472 string material_position_field_name =
"MESH_NODE_POSITIONS",
476 string element_name,
string velocity_field_name,
477 string spatial_position_field_name,
478 string material_position_field_name =
"MESH_NODE_POSITIONS",
480 Range *intersected = NULL);
483 string velocity_field_name,
string spatial_position_field_name,
484 string material_position_field_name =
"MESH_NODE_POSITIONS",
485 bool ale =
false,
bool linear =
false);
488 string velocity_field_name,
string spatial_position_field_name,
489 string material_position_field_name =
"MESH_NODE_POSITIONS",
493 string velocity_field_name,
string spatial_position_field_name,
494 string material_position_field_name =
"MESH_NODE_POSITIONS",
495 Range *forces_on_entities_ptr = NULL);
498 string velocity_field_name,
string spatial_position_field_name,
499 string material_position_field_name =
"MESH_NODE_POSITIONS",
500 bool linear =
false);
550 CHKERR MatShellGetContext(A, &void_ctx);
584 CHKERR VecAssemblyBegin(f);
593 CHKERR MatShellGetContext(A, &void_ctx);
622 CHKERR PCShellGetContext(pc, &void_ctx);
637 CHKERR PCShellGetContext(pc, &void_ctx);
675 CHKERR PCShellGetContext(pc, &void_ctx);
681 INSERT_VALUES, SCATTER_FORWARD);
683 INSERT_VALUES, SCATTER_FORWARD);
685 INSERT_VALUES, SCATTER_FORWARD);
687 INSERT_VALUES, SCATTER_FORWARD);
690 CHKERR MatMult(shell_mat_ctx->
M, shell_mat_ctx->
v,
692 CHKERR VecAXPY(shell_mat_ctx->
Ku, -1, shell_mat_ctx->
Mv);
697 CHKERR VecAXPY(shell_mat_ctx->
v, shell_mat_ctx->
ts_a,
703 INSERT_VALUES, SCATTER_REVERSE);
705 INSERT_VALUES, SCATTER_REVERSE);
707 INSERT_VALUES, SCATTER_REVERSE);
709 INSERT_VALUES, SCATTER_REVERSE);
710 CHKERR VecAssemblyBegin(x);
732#ifdef __DIRICHLET_HPP__
740 struct ShellMatrixElement :
public FEMethod {
746 typedef std::vector<PairNameFEMethodPtr> LoopsToDoType;
749 LoopsToDoType loopAuxM;
#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 ...
#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 ...
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
ublas::vector< T, ublas::bounded_array< T, N > > VectorBoundedArray
ublas::matrix< T, ublas::row_major, ublas::bounded_array< T, N > > MatrixBoundedArray
VectorBoundedArray< double, 3 > VectorDouble3
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
constexpr auto field_name
data for calculation inertia forces
Range tEts
elements in block set
double rho0
reference density
VectorDouble a0
constant acceleration
common data used by volume elements
std::vector< std::vector< double * > > jacTRowPtr
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
std::vector< std::vector< double * > > jacVelRowPtr
std::vector< VectorDouble > valVel
std::vector< MatrixDouble > jacMass
std::vector< VectorDouble > valT
std::vector< MatrixDouble > jacVel
std::vector< VectorDouble > valMass
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
std::vector< MatrixDouble > jacT
std::vector< std::vector< double * > > jacMassRowPtr
friend MoFEMErrorCode ZeroEntriesOp(Mat A)
friend MoFEMErrorCode MultOpA(Mat A, Vec x, Vec f)
definition of volume element
MoFEMErrorCode postProcess()
function is run at the end of loop
int getRule(int order)
it is used to calculate nb. of Gauss integration points
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MyVolumeFE(MoFEM::Interface &m_field)
OpEnergy(const std::string field_name, BlockData &data, CommonData &common_data, SmartPetscObj< Vec > v)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
VectorBoundedArray< adouble, 3 > a
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
MatrixBoundedArray< adouble, 9 > h
MatrixBoundedArray< adouble, 9 > g
MatrixBoundedArray< adouble, 9 > invH
MatrixBoundedArray< adouble, 9 > F
MatrixBoundedArray< adouble, 9 > H
OpEshelbyDynamicMaterialMomentumJacobian(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
MatrixBoundedArray< adouble, 9 > G
VectorBoundedArray< adouble, 3 > a_T
VectorBoundedArray< adouble, 3 > v
OpEshelbyDynamicMaterialMomentumLhs_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpEshelbyDynamicMaterialMomentumLhs_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpEshelbyDynamicMaterialMomentumLhs_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
Range forcesOnlyOnEntities
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpEshelbyDynamicMaterialMomentumRhs(const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
OpGetCommonDataAtGaussPts(const std::string field_name, CommonData &common_data)
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gardient_at_gauss_pts)
std::vector< VectorDouble > & valuesAtGaussPts
const EntityType zeroAtType
std::vector< MatrixDouble > & gradientAtGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
operator calculating deformation gradient
MatrixBoundedArray< adouble, 9 > invH
MatrixBoundedArray< adouble, 9 > G
FTensor::Index< 'k', 3 > k
FTensor::Index< 'i', 3 > i
boost::ptr_vector< MethodForForceScaling > & methodsOp
MatrixBoundedArray< adouble, 9 > H
OpMassJacobian(const std::string field_name, BlockData &data, CommonData &common_data, boost::ptr_vector< MethodForForceScaling > &methods_op, int tag, bool linear=false)
MatrixBoundedArray< adouble, 9 > g
MatrixBoundedArray< adouble, 9 > h
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
VectorBoundedArray< adouble, 3 > a_res
VectorBoundedArray< adouble, 3 > dp_dt
MatrixBoundedArray< adouble, 9 > F
VectorBoundedArray< adouble, 3 > dot_W
std::vector< double > active
FTensor::Index< 'j', 3 > j
VectorBoundedArray< adouble, 3 > a
OpMassLhs_dM_dX(const std::string field_name, const std::string col_field, BlockData &data, CommonData &common_data)
MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpMassLhs_dM_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr=NULL)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
Range forcesOnlyOnEntities
MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpMassLhs_dM_dx(const std::string field_name, const std::string col_field, BlockData &data, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpMassRhs(const std::string field_name, BlockData &data, CommonData &common_data)
VectorBoundedArray< adouble, 3 > v
VectorBoundedArray< adouble, 3 > dot_u
VectorBoundedArray< adouble, 3 > a_res
std::vector< double > active
MatrixBoundedArray< adouble, 9 > H
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpVelocityJacobian(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
VectorBoundedArray< adouble, 3 > dot_w
VectorBoundedArray< adouble, 3 > dot_W
MatrixBoundedArray< adouble, 9 > h
MatrixBoundedArray< adouble, 9 > F
MatrixBoundedArray< adouble, 9 > invH
OpVelocityLhs_dV_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpVelocityLhs_dV_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpVelocityLhs_dV_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpVelocityRhs(const std::string field_name, BlockData &data, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
friend MoFEMErrorCode PCShellSetUpOp(PC pc)
friend MoFEMErrorCode PCShellDestroy(PC pc)
bool initPC
check if PC is initialized
PCShellCtx(Mat shell_mat)
friend MoFEMErrorCode PCShellApplyOp(PC pc, Vec f, Vec x)
MatShellCtx * shellMatCtx
pointer to shell matrix
MoFEM::Interface & mField
MoFEMErrorCode preProcess()
Calculate inconsistency between approximation of velocities and velocities calculated from displaceme...
ShellResidualElement(MoFEM::Interface &m_field)
const std::string spatialPositionField
const std::string velocityField
UpdateAndControl(MoFEM::Interface &m_field, TS _ts, const std::string velocity_field, const std::string spatial_position_field)
MoFEM::Interface & mField
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode preProcess()
Scatter values from t_u_dt on the fields.
structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_e...
static MoFEMErrorCode MultOpA(Mat A, Vec x, Vec f)
Mult operator for shell matrix.
ConvectiveMassElement(MoFEM::Interface &m_field, short int tag)
static MoFEMErrorCode ZeroEntriesOp(Mat A)
MoFEMErrorCode setVelocityOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false)
MyVolumeFE feEnergy
calculate kinetic energy
MoFEMErrorCode setShellMatrixMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool linear=false)
MyVolumeFE feVelRhs
calculate right hand side for tetrahedral elements
MoFEMErrorCode addVelocityElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
MyVolumeFE & getLoopFeVelRhs()
get rhs volume element
static MoFEMErrorCode PCShellDestroy(PC pc)
MoFEMErrorCode addConvectiveMassElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
MyVolumeFE & getLoopFeTRhs()
get rhs volume element
static MoFEMErrorCode PCShellApplyOp(PC pc, Vec f, Vec x)
apply pre-conditioner for shell matrix
MoFEMErrorCode setKinematicEshelbyOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", Range *forces_on_entities_ptr=NULL)
MoFEMErrorCode addEshelbyDynamicMaterialMomentum(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel(), Range *intersected=NULL)
MoFEMErrorCode addHOOpsVol()
MyVolumeFE feTRhs
calculate right hand side for tetrahedral elements
MyVolumeFE feMassRhs
calculate right hand side for tetrahedral elements
MyVolumeFE & getLoopFeMassRhs()
get rhs volume element
MyVolumeFE feTLhs
calculate left hand side for tetrahedral elements
MyVolumeFE & getLoopFeEnergy()
get kinetic energy element
MoFEMErrorCode setConvectiveMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, bool linear=false)
boost::ptr_vector< MethodForForceScaling > methodsOp
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
MyVolumeFE & getLoopFeMassLhs()
get lhs volume element
static MoFEMErrorCode PCShellSetUpOp(PC pc)
MoFEM::Interface & mField
MyVolumeFE & getLoopFeVelLhs()
get lhs volume element
MoFEMErrorCode setBlocks()
MyVolumeFE feVelLhs
calculate left hand side for tetrahedral elements
MyVolumeFE & getLoopFeMassAuxLhs()
get lhs volume element for Kuu shell matrix
MyVolumeFE & getLoopFeTLhs()
get lhs volume element
Set Dirichlet boundary conditions on displacements.
Deprecated interface functions.
structure for User Loop Methods on finite elements
intrusive_ptr for managing petsc objects
Volume finite element base.