9#ifndef __CONVECTIVE_MASS_ELEMENT_HPP
10#define __CONVECTIVE_MASS_ELEMENT_HPP
13#error "MoFEM need to be compiled with ADOL-C"
31 struct MyVolumeFE :
public VolumeElementForcesAndSourcesCore {
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);
236 virtual MoFEMErrorCode
getJac(EntitiesFieldData::EntData &col_data,
241 EntitiesFieldData::EntData &row_data,
242 EntitiesFieldData::EntData &col_data);
250 MoFEMErrorCode
getJac(EntitiesFieldData::EntData &col_data,
int gg);
258 MoFEMErrorCode
getJac(EntitiesFieldData::EntData &col_data,
int gg);
261 struct OpEnergy :
public VolumeElementForcesAndSourcesCore::UserDataOperator,
266 SmartPetscObj<Vec>
V;
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);
323 virtual MoFEMErrorCode
getJac(EntitiesFieldData::EntData &col_data,
333 virtual MoFEMErrorCode
getJac(EntitiesFieldData::EntData &col_data,
343 virtual MoFEMErrorCode
getJac(EntitiesFieldData::EntData &col_data,
348 :
public VolumeElementForcesAndSourcesCore::UserDataOperator,
360 bool jacobian =
true);
362 VectorBoundedArray<adouble, 3>
a,
v,
a_T;
367 EntitiesFieldData::EntData &row_data);
371 :
public VolumeElementForcesAndSourcesCore::UserDataOperator,
381 Range *forcesonlyonentities_ptr);
386 EntitiesFieldData::EntData &row_data);
395 Range *forcesonlyonentities_ptr);
397 virtual MoFEMErrorCode
getJac(EntitiesFieldData::EntData &col_data,
408 Range *forcesonlyonentities_ptr);
410 virtual MoFEMErrorCode
getJac(EntitiesFieldData::EntData &col_data,
421 Range *forcesonlyonentities_ptr);
423 virtual MoFEMErrorCode
getJac(EntitiesFieldData::EntData &col_data,
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);
514 MoFEMErrorCode
iNit();
518 friend MoFEMErrorCode
MultOpA(Mat
A, Vec x, Vec f);
550 CHKERR MatShellGetContext(
A, &void_ctx);
584 CHKERR VecAssemblyBegin(f);
593 CHKERR MatShellGetContext(
A, &void_ctx);
609 MoFEMErrorCode
iNit();
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.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
UBlasVector< double > VectorDouble
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()
int getRule(int order)
it is used to calculate nb. of Gauss integration points
MoFEMErrorCode preProcess()
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
MatrixBoundedArray< adouble, 9 > G
VectorBoundedArray< adouble, 3 > a_T
VectorBoundedArray< adouble, 3 > v
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
Range forcesOnlyOnEntities
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
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
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
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)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
Range forcesOnlyOnEntities
MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_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)
VectorBoundedArray< adouble, 3 > dot_w
VectorBoundedArray< adouble, 3 > dot_W
MatrixBoundedArray< adouble, 9 > h
MatrixBoundedArray< adouble, 9 > F
MatrixBoundedArray< adouble, 9 > invH
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
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...
const std::string spatialPositionField
const std::string velocityField
MoFEM::Interface & mField
MoFEMErrorCode postProcess()
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.
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.