10using namespace boost::numeric;
13static char help[] =
"-order approximation order\n"
16struct OpK :
public VolumeElementForcesAndSourcesCore::UserDataOperator {
64 tD(
i,
j,
k,
l) *= coefficient;
106 if (row_side == col_side && row_type == col_type) {
139 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
141 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
142 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
143 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2));
152 double vol = getVolume();
155 auto t_w = getFTensor0IntegrationWeight();
163 const double a = t_w * vol;
166 for (
int rr = 0; rr !=
nbRows / 3; ++rr) {
169 auto t_m = get_tensor2(
K, 3 * rr, 0);
175 for (
int cc = 0; cc !=
nbCols / 3; ++cc) {
179 a * (
tD(
i,
j,
k,
l) * (t_row_diff_base(
j) * t_col_diff_base(
l)));
209 const int *row_indices = &*row_data.
getIndices().data().begin();
211 const int *col_indices = &*col_data.
getIndices().data().begin();
212 Mat
B = getFEMethod()->ksp_B != PETSC_NULLPTR ? getFEMethod()->ksp_B
213 : getFEMethod()->snes_B;
216 &*
K.data().begin(), ADD_VALUES);
223 &*
K.data().begin(), ADD_VALUES);
254 nF.resize(nb_dofs,
false);
258 const int nb_gauss_pts = data.
getN().size1();
271 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
273 double w = 0.5 * t_w;
279 for (
int bb = 0; bb != nb_dofs / 3; ++bb) {
306 const Range &fix_second_node)
315 std::set<int> set_fix_dofs;
318 if (dit->get()->getDofCoeffIdx() == 2) {
320 set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());
325 if (dit->get()->getDofCoeffIdx() == 1) {
326 set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());
331 set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());
335 std::vector<int> fix_dofs(set_fix_dofs.size());
337 std::copy(set_fix_dofs.begin(), set_fix_dofs.end(), fix_dofs.begin());
339 CHKERR MatAssemblyBegin(
ksp_B, MAT_FINAL_ASSEMBLY);
348 CHKERR MatZeroRowsColumns(
ksp_B, fix_dofs.size(), &fix_dofs[0], 1,
x,
369 int operator()(
int,
int,
int p)
const {
return 2 * (p - 1); }
390int main(
int argc,
char *argv[]) {
392 const string default_options =
"-ksp_type fgmres \n"
394 "-pc_factor_mat_solver_type mumps \n"
395 "-mat_mumps_icntl_20 0 \n"
398 string param_file =
"param_file.petsc";
399 if (!
static_cast<bool>(ifstream(param_file))) {
400 std::ofstream file(param_file.c_str(), std::ios::ate);
401 if (file.is_open()) {
402 file << default_options;
410 moab::Core mb_instance;
411 moab::Interface &moab = mb_instance;
422 PetscBool flg_test = PETSC_FALSE;
423 PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"SimpleElasticProblem",
431 CHKERR PetscOptionsBool(
"-test",
"if true is ctest",
"", flg_test,
432 &flg_test, PETSC_NULLPTR);
441 Range fix_faces, pressure_faces, fix_nodes, fix_second_node;
446 BRICK_PRESSURE_FACES = 3,
452 int id =
bit->getMeshsetId();
454 if (
id == FIX_BRICK_FACES) {
460 CHKERR m_field.
get_moab().get_adjacencies(fix_faces, 0,
false, adj_ents,
461 moab::Interface::UNION);
463 CHKERR m_field.
get_moab().get_adjacencies(fix_faces, 1,
false, adj_ents,
464 moab::Interface::UNION);
465 fix_faces.merge(adj_ents);
466 }
else if (
id == FIX_NODES) {
471 }
else if (
id == BRICK_PRESSURE_FACES) {
473 meshset, 2, pressure_faces,
true);
478 meshset, 0, fix_second_node,
true);
515 boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(
517 boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(
521 elastic_fe->getRuleHook =
VolRule();
522 pressure_fe->getRuleHook =
FaceRule();
526 elastic_fe->getOpPtrVector().push_back(
new OpK());
528 pressure_fe->getOpPtrVector().push_back(
new OpPressure());
530 boost::shared_ptr<FEMethod> fix_dofs_fe(
533 boost::shared_ptr<FEMethod> null_fe;
537 dm, simple_interface->
getDomainFEName(), elastic_fe, null_fe, null_fe);
549 CHKERR DMCreateMatrix(dm, &A);
552 CHKERR DMCreateGlobalVector(dm, &x);
556 CHKERR VecDuplicate(x, &f);
560 elastic_fe->ksp_B =
A;
561 fix_dofs_fe->ksp_B =
A;
565 fix_dofs_fe->ksp_f = f;
566 pressure_fe->ksp_f = f;
582 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
586 CHKERR KSPSetFromOptions(solver);
589 CHKERR KSPSetOperators(solver, A, A);
595 CHKERR KSPSolve(solver, f, x);
598 CHKERR KSPDestroy(&solver);
602 VecView(x, PETSC_VIEWER_STDOUT_WORLD);
609 auto get_post_proc_ele = [&]() {
610 auto jac_ptr = boost::make_shared<MatrixDouble>();
611 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
612 auto det_ptr = boost::make_shared<VectorDouble>();
614 auto post_proc_ele = boost::make_shared<
618 post_proc_ele->getOpPtrVector().push_back(
620 post_proc_ele->getOpPtrVector().push_back(
622 post_proc_ele->getOpPtrVector().push_back(
625 auto u_ptr = boost::make_shared<MatrixDouble>();
626 auto grad_ptr = boost::make_shared<MatrixDouble>();
627 post_proc_ele->getOpPtrVector().push_back(
629 post_proc_ele->getOpPtrVector().push_back(
634 post_proc_ele->getOpPtrVector().push_back(
638 post_proc_ele->getPostProcMesh(), post_proc_ele->getMapGaussPts(),
643 {{
"GRAD", grad_ptr}},
649 return post_proc_ele;
652 if (
auto post_proc = get_post_proc_ele()) {
656 CHKERR post_proc->writeFile(
"out.h5m");
660 if (flg_test == PETSC_TRUE) {
661 const double x_vec_norm_const = 0.4;
666 CHKERR VecNorm(x, NORM_INFINITY, &norm_check);
667 if (std::abs(norm_check - x_vec_norm_const) / x_vec_norm_const >
670 "test failed (nrm 2 %6.4e)", norm_check);
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#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_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set compute operator for KSP solver via sub-matrix and IS.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set KSP operators and push mofem finite element methods.
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_NUMEREDDOF_ROW_FOR_LOOP_(PROBLEMPTR, IT)
use with loops to iterate row DOFs
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FTensor::Index< 'm', 3 > m
ApplyDirichletBc(const Range &fix_faces, const Range &fix_nodes, const Range &fix_second_node)
MoFEMErrorCode postProcess()
function is run at the end of loop
Set integration rule to boundary elements.
int operator()(int, int, int p) const
const Problem * problemPtr
raw pointer to problem
virtual moab::Interface & get_moab()=0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
structure for User Loop Methods on finite elements
default operator for TRI element
auto getFTensor1Normal()
get normal as tensor
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Get values at integration pts for tensor field rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
Simple interface for fast problem set-up.
MoFEMErrorCode buildProblem()
Build problem.
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
MoFEMErrorCode defineFiniteElements()
Define finite elements.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
MoFEMErrorCode buildFiniteElements()
Build finite elements.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode buildFields()
Build fields.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
const std::string getDomainFEName() const
Get the Domain FE Name.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Volume finite element base.
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate B^T D B operator.
bool isDiag
true if this block is on diagonal
FTensor::Ddg< double, 3, 3 > tD
int nbIntegrationPts
number of integration points
int nbRows
number of dofs on rows
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations for give operator.
int nbCols
number if dof on column
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble local entity block matrix.
FTensor::Index< 'i', 3 > i
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpPressure(const double pressure_val=1)
Set integration rule to volume elements.
int operator()(int, int, int p) const