21 #include <BasicFiniteElements.hpp>
22 using namespace MoFEM;
24 static char help[] =
"...\n\n";
32 const double D = 2e-3;
36 const double u0 = 0.1;
62 :
OpEle(
"u",
"u",
OpEle::OPROWCOL), commonData(data) {
66 EntityType col_type,
EntData &row_data,
69 const int nb_row_dofs = row_data.
getIndices().size();
70 const int nb_col_dofs = col_data.
getIndices().size();
71 if (nb_row_dofs && nb_col_dofs) {
72 const int nb_integration_pts = getGaussPts().size2();
73 mat.resize(nb_row_dofs, nb_col_dofs,
false);
76 auto t_w = getFTensor0IntegrationWeight();
77 const double vol = getMeasure();
78 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
79 const double a = t_w * vol;
80 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
82 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
83 mat(rr, cc) += a * t_row_base * t_col_base;
93 if (row_side != col_side || row_type != col_type) {
94 transMat.resize(nb_col_dofs, nb_row_dofs,
false);
95 noalias(transMat) = trans(mat);
117 :
OpEle(
"u",
OpEle::OPROW), commonData(data) {}
122 vecF.resize(nb_dofs,
false);
125 const int nb_integration_pts = getGaussPts().size2();
128 auto t_w = getFTensor0IntegrationWeight();
130 const double vol = getMeasure();
131 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
132 const double a = vol * t_w;
133 const double f = a *
r * t_val * (1 - t_val /
k);
134 for (
int rr = 0; rr != nb_dofs; ++rr) {
135 const double b =
f * t_base;
144 CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
166 :
OpEle(
"u",
OpEle::OPROW), commonData(data) {}
171 vecF.resize(nb_dofs,
false);
174 const int nb_integration_pts = getGaussPts().size2();
176 auto t_grad = getFTensor1FromMat<DIM>(commonData->grad);
179 auto t_w = getFTensor0IntegrationWeight();
182 const double vol = getMeasure();
183 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
184 const double a = vol * t_w;
185 for (
int rr = 0; rr != nb_dofs; ++rr) {
186 vecF[rr] += a * (t_base * t_dot_val +
D * t_diff_base(
i) * t_grad(
i));
195 CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
218 :
OpEle(
"u",
"u",
OpEle::OPROWCOL), commonData(data) {
222 EntityType col_type,
EntData &row_data,
225 const int nb_row_dofs = row_data.
getIndices().size();
226 const int nb_col_dofs = col_data.
getIndices().size();
227 if (nb_row_dofs && nb_col_dofs) {
228 mat.resize(nb_row_dofs, nb_col_dofs,
false);
231 const int nb_integration_pts = getGaussPts().size2();
234 auto t_w = getFTensor0IntegrationWeight();
237 const double ts_a = getFEMethod()->ts_a;
238 const double vol = getMeasure();
239 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
240 const double a = vol * t_w;
241 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
246 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
247 mat(rr, cc) += a * (t_row_base * t_col_base * ts_a +
248 D * t_row_diff_base(
i) * t_col_diff_base(
i));
261 if (row_side != col_side || row_type != col_type) {
262 transMat.resize(nb_col_dofs, nb_row_dofs,
false);
263 noalias(transMat) = trans(mat);
265 &transMat(0, 0), ADD_VALUES);
285 boost::shared_ptr<PostProcFaceOnRefinedMesh> &post_proc)
286 : dM(dm), postProc(post_proc){};
295 CHKERR postProc->writeFile(
296 "out_level_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
303 boost::shared_ptr<PostProcFaceOnRefinedMesh>
postProc;
310 int main(
int argc,
char *argv[]) {
325 DMType dm_name =
"DMMOFEM";
346 boost::shared_ptr<CommonData> data(
new CommonData());
348 auto val_ptr = boost::shared_ptr<VectorDouble>(data, &data->val);
349 auto dot_val_ptr = boost::shared_ptr<VectorDouble>(data, &data->dot_val);
350 auto grad_ptr = boost::shared_ptr<MatrixDouble>(data, &data->grad);
354 boost::shared_ptr<Ele> vol_ele_slow_rhs(
new Ele(m_field));
355 boost::shared_ptr<Ele> vol_ele_stiff_rhs(
new Ele(m_field));
356 boost::shared_ptr<Ele> vol_ele_stiff_lhs(
new Ele(m_field));
359 vol_ele_slow_rhs->getOpPtrVector().push_back(
366 auto solve_for_g = [&]() {
368 if (vol_ele_slow_rhs->vecAssembleSwitch) {
369 CHKERR VecGhostUpdateBegin(vol_ele_slow_rhs->ts_F, ADD_VALUES,
371 CHKERR VecGhostUpdateEnd(vol_ele_slow_rhs->ts_F, ADD_VALUES,
373 CHKERR VecAssemblyBegin(vol_ele_slow_rhs->ts_F);
374 CHKERR VecAssemblyEnd(vol_ele_slow_rhs->ts_F);
375 *vol_ele_slow_rhs->vecAssembleSwitch =
false;
377 CHKERR KSPSolve(data->ksp, vol_ele_slow_rhs->ts_F,
378 vol_ele_slow_rhs->ts_F);
382 vol_ele_slow_rhs->postProcessHook = solve_for_g;
385 vol_ele_stiff_rhs->getOpPtrVector().push_back(
387 vol_ele_stiff_rhs->getOpPtrVector().push_back(
389 vol_ele_stiff_rhs->getOpPtrVector().push_back(
391 vol_ele_stiff_rhs->getOpPtrVector().push_back(
393 vol_ele_stiff_rhs->getOpPtrVector().push_back(
397 vol_ele_stiff_lhs->getOpPtrVector().push_back(
399 vol_ele_stiff_lhs->getOpPtrVector().push_back(
401 vol_ele_stiff_lhs->getOpPtrVector().push_back(
405 auto vol_rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
406 vol_ele_slow_rhs->getRuleHook = vol_rule;
407 vol_ele_stiff_rhs->getRuleHook = vol_rule;
408 vol_ele_stiff_lhs->getRuleHook = vol_rule;
411 boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc =
412 boost::shared_ptr<PostProcFaceOnRefinedMesh>(
414 boost::shared_ptr<ForcesAndSourcesCore>
null;
416 post_proc->generateReferenceElementMesh();
418 post_proc->addFieldValuesPostProc(
"u");
421 auto dm = simple_interface->
getDM();
429 1,
BLOCKSET, 2, inner_surface,
true);
430 if (!inner_surface.empty()) {
431 Range inner_surface_verts;
432 CHKERR moab.get_connectivity(inner_surface, inner_surface_verts,
false);
434 u0, MBVERTEX, inner_surface_verts,
"u");
441 CHKERR moab.get_entities_by_dimension(0, 2, surface,
false);
444 CHKERR skin.find_skin(0, surface,
false, edges);
446 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
447 CHKERR pcomm->filter_pstatus(edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
448 PSTATUS_NOT, -1, &edges_part);
450 CHKERR moab.get_connectivity(edges_part, edges_verts,
false);
455 unite(edges_verts, edges_part));
459 CHKERR MatZeroEntries(data->M);
460 boost::shared_ptr<Ele> vol_mass_ele(
new Ele(m_field));
461 vol_mass_ele->getOpPtrVector().push_back(
new OpAssembleMass(data));
464 CHKERR MatAssemblyBegin(data->M, MAT_FINAL_ASSEMBLY);
465 CHKERR MatAssemblyEnd(data->M, MAT_FINAL_ASSEMBLY);
470 CHKERR KSPSetOperators(data->ksp, data->M, data->M);
471 CHKERR KSPSetFromOptions(data->ksp);
472 CHKERR KSPSetUp(data->ksp);
477 CHKERR TSSetType(ts, TSARKIMEX);
478 CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
482 vol_ele_stiff_lhs,
null,
null);
485 vol_ele_stiff_rhs,
null,
null);
488 vol_ele_slow_rhs,
null,
null);
491 boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(dm, post_proc));
493 monitor_ptr,
null,
null);
503 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
504 CHKERR TSSetSolution(ts, X);
505 CHKERR TSSetFromOptions(ts);
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MYPCOMM_INDEX
default communicator number PCOMM
#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.
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
ublas::vector< double, DoubleAllocator > VectorDouble
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
implementation of Data Operators for Forces and Sources
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
PetscErrorCode DMMoFEMTSSetRHSFunction(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS the right hand side function
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
OpCalculateInvJacForFaceImpl< 2 > OpCalculateInvJacForFace
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
DeprecatedCoreInterface Interface
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
const double D
diffusivity
DataForcesAndSourcesCore::EntData EntData
const int save_every_nth_step
const double r
rate factor
const double u0
inital vale on blocksets
FaceElementForcesAndSourcesCore Ele
int main(int argc, char *argv[])
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =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.
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
FTensor::Tensor1< double *, 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.
Deprecated interface functions.
structure for User Loop Methods on finite elements
default operator for TRI element
Interface for managing meshsets containing materials and boundary conditions.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get rate of scalar field at integration points.
Get value at integration points for scalar field.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
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 getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getProblemName() const
Get the Problem Name.
const std::string getDomainFEName() const
Get the Domain FE Name.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
SmartPetscObj< KSP > ksp
Linear solver.
VectorDouble dot_val
Rate of values of field "u" at integration points.
VectorDouble val
Values of field "u" at integration points.
SmartPetscObj< Mat > M
Mass matrix.
MatrixDouble grad
Gradients of field "u" at integration points.
MatrixDouble invJac
Inverse of element jacobian.
boost::shared_ptr< PostProcFaceOnRefinedMesh > postProc
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode operator()()
function is run for every finite element
Monitor(SmartPetscObj< DM > &dm, boost::shared_ptr< PostProcFaceOnRefinedMesh > &post_proc)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
boost::shared_ptr< CommonData > commonData
OpAssembleMass(boost::shared_ptr< CommonData > &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpAssembleSlowRhs(boost::shared_ptr< CommonData > &data)
boost::shared_ptr< CommonData > commonData
Assemble stiff part tangent.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
boost::shared_ptr< CommonData > commonData
OpAssembleStiffLhs(boost::shared_ptr< CommonData > &data)
OpAssembleStiffRhs(boost::shared_ptr< CommonData > &data)
boost::shared_ptr< CommonData > commonData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.