10static char help[] =
"...\n\n";
56 const int nb_row_dofs = row_data.
getIndices().size();
57 const int nb_col_dofs = col_data.
getIndices().size();
58 if (nb_row_dofs && nb_col_dofs) {
59 const int nb_integration_pts =
getGaussPts().size2();
60 mat.resize(nb_row_dofs, nb_col_dofs,
false);
65 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
66 const double a = t_w * vol;
67 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
69 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
70 mat(rr, cc) +=
a * t_row_base * t_col_base;
80 if (row_side != col_side || row_type != col_type) {
81 transMat.resize(nb_col_dofs, nb_row_dofs,
false);
109 vecF.resize(nb_dofs,
false);
112 const int nb_integration_pts =
getGaussPts().size2();
118 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
119 const double a = vol * t_w;
120 const double f =
a *
r * t_val * (1 - t_val /
k);
121 for (
int rr = 0; rr != nb_dofs; ++rr) {
122 const double b = f * t_base;
158 vecF.resize(nb_dofs,
false);
161 const int nb_integration_pts =
getGaussPts().size2();
163 auto t_grad = getFTensor1FromMat<DIM>(
commonData->grad);
170 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
171 const double a = vol * t_w;
172 for (
int rr = 0; rr != nb_dofs; ++rr) {
173 vecF[rr] +=
a * (t_base * t_dot_val +
D * t_diff_base(
i) * t_grad(
i));
212 const int nb_row_dofs = row_data.
getIndices().size();
213 const int nb_col_dofs = col_data.
getIndices().size();
214 if (nb_row_dofs && nb_col_dofs) {
215 mat.resize(nb_row_dofs, nb_col_dofs,
false);
218 const int nb_integration_pts =
getGaussPts().size2();
226 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
227 const double a = vol * t_w;
228 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
233 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
234 mat(rr, cc) +=
a * (t_row_base * t_col_base * ts_a +
235 D * t_row_diff_base(
i) * t_col_diff_base(
i));
248 if (row_side != col_side || row_type != col_type) {
249 transMat.resize(nb_col_dofs, nb_row_dofs,
false);
282 "out_level_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
296int main(
int argc,
char *argv[]) {
305 moab::Core mb_instance;
306 moab::Interface &moab = mb_instance;
311 DMType dm_name =
"DMMOFEM";
332 boost::shared_ptr<CommonData> data(
new CommonData());
334 auto val_ptr = boost::shared_ptr<VectorDouble>(data, &data->val);
335 auto dot_val_ptr = boost::shared_ptr<VectorDouble>(data, &data->dot_val);
336 auto grad_ptr = boost::shared_ptr<MatrixDouble>(data, &data->grad);
340 boost::shared_ptr<Ele> vol_ele_slow_rhs(
new Ele(m_field));
341 boost::shared_ptr<Ele> vol_ele_stiff_rhs(
new Ele(m_field));
342 boost::shared_ptr<Ele> vol_ele_stiff_lhs(
new Ele(m_field));
345 vol_ele_slow_rhs->getOpPtrVector().push_back(
352 auto solve_for_g = [&]() {
354 if (*(vol_ele_slow_rhs->vecAssembleSwitch)) {
355 CHKERR VecGhostUpdateBegin(vol_ele_slow_rhs->ts_F, ADD_VALUES,
357 CHKERR VecGhostUpdateEnd(vol_ele_slow_rhs->ts_F, ADD_VALUES,
359 CHKERR VecAssemblyBegin(vol_ele_slow_rhs->ts_F);
360 CHKERR VecAssemblyEnd(vol_ele_slow_rhs->ts_F);
361 *vol_ele_slow_rhs->vecAssembleSwitch =
false;
363 CHKERR KSPSolve(data->ksp, vol_ele_slow_rhs->ts_F,
364 vol_ele_slow_rhs->ts_F);
368 vol_ele_slow_rhs->postProcessHook = solve_for_g;
370 auto det_ptr = boost::make_shared<VectorDouble>();
371 auto jac_ptr = boost::make_shared<MatrixDouble>();
372 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
374 vol_ele_stiff_rhs->getOpPtrVector().push_back(
376 vol_ele_stiff_rhs->getOpPtrVector().push_back(
378 vol_ele_stiff_rhs->getOpPtrVector().push_back(
380 vol_ele_stiff_rhs->getOpPtrVector().push_back(
382 vol_ele_stiff_rhs->getOpPtrVector().push_back(
384 vol_ele_stiff_rhs->getOpPtrVector().push_back(
388 vol_ele_stiff_lhs->getOpPtrVector().push_back(
390 vol_ele_stiff_lhs->getOpPtrVector().push_back(
392 vol_ele_stiff_lhs->getOpPtrVector().push_back(
394 vol_ele_stiff_lhs->getOpPtrVector().push_back(
398 auto vol_rule = [](int, int,
int p) ->
int {
return 2 *
p; };
399 vol_ele_slow_rhs->getRuleHook = vol_rule;
400 vol_ele_stiff_rhs->getRuleHook = vol_rule;
401 vol_ele_stiff_lhs->getRuleHook = vol_rule;
405 auto post_proc = boost::make_shared<PostProcEle>(m_field);
406 boost::shared_ptr<ForcesAndSourcesCore> null;
410 auto u_ptr = boost::make_shared<VectorDouble>();
411 post_proc->getOpPtrVector().push_back(
413 post_proc->getOpPtrVector().push_back(
417 post_proc->getPostProcMesh(), post_proc->getMapGaussPts(),
432 auto dm = simple_interface->getDM();
440 1,
BLOCKSET, 2, inner_surface,
true);
441 if (!inner_surface.empty()) {
442 Range inner_surface_verts;
443 CHKERR moab.get_connectivity(inner_surface, inner_surface_verts,
false);
445 u0, MBVERTEX, inner_surface_verts,
"u");
453 Skinner skin(&m_field.get_moab());
457 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
458 CHKERR pcomm->filter_pstatus(edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
459 PSTATUS_NOT, -1, &edges_part);
461 CHKERR moab.get_connectivity(edges_part, edges_verts,
false);
465 simple_interface->getProblemName(),
"u",
466 unite(edges_verts, edges_part));
470 CHKERR MatZeroEntries(data->M);
471 boost::shared_ptr<Ele> vol_mass_ele(
new Ele(m_field));
472 vol_mass_ele->getOpPtrVector().push_back(
new OpAssembleMass(data));
475 CHKERR MatAssemblyBegin(data->M, MAT_FINAL_ASSEMBLY);
476 CHKERR MatAssemblyEnd(data->M, MAT_FINAL_ASSEMBLY);
480 data->ksp =
createKSP(m_field.get_comm());
481 CHKERR KSPSetOperators(data->ksp, data->M, data->M);
482 CHKERR KSPSetFromOptions(data->ksp);
483 CHKERR KSPSetUp(data->ksp);
486 auto ts =
createTS(m_field.get_comm());
488 CHKERR TSSetType(ts, TSARKIMEX);
489 CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
493 vol_ele_stiff_lhs, null, null);
496 vol_ele_stiff_rhs, null, null);
499 vol_ele_slow_rhs, null, null);
502 boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(dm, post_proc));
504 monitor_ptr, null, null);
514 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
515 CHKERR TSSetSolution(ts, X);
516 CHKERR TSSetFromOptions(ts);
#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 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 DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
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.
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
auto createKSP(MPI_Comm comm)
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)
auto createTS(MPI_Comm comm)
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
const double D
diffusivity
const int save_every_nth_step
const double r
rate factor
const double u0
inital vale on blocksets
const double k
caring capacity
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FaceElementForcesAndSourcesCore Ele
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.
bool sYmm
If true assume that matrix is symmetric structure.
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.
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 getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume 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 value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
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 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.
intrusive_ptr for managing petsc objects
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
[Push operators to pipeline]
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.
boost::shared_ptr< PostProcEle > 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< PostProcEle > &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.