8static char help[] =
"...\n\n";
24 double operator()(
const double x,
const double y,
const double z,
25 const double t)
const {
26 if (y <= -1.6000 &&
t <= 0.50000) {
36 return factor * (1.0 * (
c * u * (u -
alpha) * (1.0 - u) - u *
v));
59const double Tol = 1e-9;
64 double operator()(
const double u,
const double v,
const double dt)
const {
74 double rhsv =
rhs_v(u, vk);
75 double Drhsv =
Drhs_v(u, vk);
76 dv = - 1.0 / (1 -
dt * Drhsv) * (vk -
v -
dt * rhsv);
101 post_proc = boost::shared_ptr<PostProcFaceOnRefinedMesh>(
109 boost::shared_ptr<MatrixDouble>(
data1, &
data1->flux_values);
114 boost::shared_ptr<VectorDouble>(
data1, &
data1->mass_values);
119 boost::shared_ptr<MatrixDouble>(
data2, &
data2->flux_values);
123 boost::shared_ptr<VectorDouble>(
data2, &
data2->mass_values);
141 boost::shared_ptr<VectorDouble> &mass_ptr);
143 boost::shared_ptr<PreviousData> &
data1,
144 boost::shared_ptr<PreviousData> &
data2);
146 boost::shared_ptr<PreviousData> &
data);
149 boost::shared_ptr<VectorDouble> &mass_ptr,
150 boost::shared_ptr<MatrixDouble> &flux_ptr,
151 boost::shared_ptr<VectorDouble> &mass_dot_ptr,
152 boost::shared_ptr<VectorDouble> &flux_div_ptr);
154 boost::shared_ptr<PreviousData> &
data1,
155 boost::shared_ptr<PreviousData> &
data2,
156 std::map<int, BlockData> &block_map,
159 std::string flux_field,
160 boost::shared_ptr<VectorDouble> &mass_ptr,
161 boost::shared_ptr<MatrixDouble> &flux_ptr);
163 boost::shared_ptr<PreviousData> &datau,
164 boost::shared_ptr<PreviousData> &datav,
165 std::map<int, BlockData> &block_map);
169 boost::shared_ptr<FaceEle> &initial_ele,
205 boost::shared_ptr<PreviousData>
data1;
206 boost::shared_ptr<PreviousData>
data2;
225 boost::shared_ptr<ForcesAndSourcesCore>
null;
262 string name = it->getName();
263 const int id = it->getMeshsetId();
264 if (name.compare(0, 14,
"REGION1") == 0) {
266 id,
BLOCKSET, 2, block_map[
id].block_ents,
true);
267 block_map[id].B0 = 1e-3;
269 }
else if (name.compare(0, 14,
"REGION2") == 0) {
271 id,
BLOCKSET, 2, block_map[
id].block_ents,
true);
272 block_map[id].B0 = 1e-1;
280 std::string natural) {
283 string name = it->getName();
284 if (name.compare(0, 14, natural) == 0) {
288 }
else if (name.compare(0, 14, essential) == 0) {
308 boost::shared_ptr<VectorDouble> &mass_ptr) {
317 boost::shared_ptr<PreviousData> &data_1,
318 boost::shared_ptr<PreviousData> &data_2) {
334 boost::shared_ptr<PreviousData> &data) {
336 auto det_ptr = boost::make_shared<VectorDouble>();
337 auto jac_ptr = boost::make_shared<MatrixDouble>();
338 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
340 vol_ele->getOpPtrVector().push_back(
343 vol_ele->getOpPtrVector().push_back(
351 boost::shared_ptr<VectorDouble> &mass_ptr,
352 boost::shared_ptr<MatrixDouble> &flux_ptr,
353 boost::shared_ptr<VectorDouble> &mass_dot_ptr,
354 boost::shared_ptr<VectorDouble> &flux_div_ptr) {
373 std::string flux_field,
374 boost::shared_ptr<PreviousData> &data1,
375 boost::shared_ptr<PreviousData> &data2,
376 std::map<int, BlockData> &block_map,
389 boost::shared_ptr<VectorDouble> &mass_ptr,
390 boost::shared_ptr<MatrixDouble> &flux_ptr) {
401 std::string flux_field,
402 boost::shared_ptr<PreviousData> &data1,
403 boost::shared_ptr<PreviousData> &data2,
404 std::map<int, BlockData> &block_map) {
423 auto vol_rule = [](int, int,
int p) ->
int {
return 2 *
p; };
434 boost::shared_ptr<FaceEle> &initial_ele,
437 initial_ele->getOpPtrVector().push_back(
452 CHKERR TSARKIMEXSetType(
ts, TSARKIMEXA2);
491 CHKERR TSSetDuration(
ts, PETSC_DEFAULT, ftime);
499 CHKERR SNESGetKSP(snes, &ksp);
501 CHKERR KSPGetPC(ksp, &pc);
502 PetscBool is_pcfs = PETSC_FALSE;
503 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
506 if (is_pcfs == PETSC_TRUE) {
511 problem_ptr->
getName(),
ROW,
"U", 0, 1, &is_mass);
513 problem_ptr->
getName(),
ROW,
"F", 0, 1, &is_flux);
517 CHKERR PCFieldSplitSetIS(pc, NULL, is_mass);
518 CHKERR PCFieldSplitSetIS(pc, NULL, is_flux);
521 CHKERR ISDestroy(&is_flux);
522 CHKERR ISDestroy(&is_mass);
579 boost::shared_ptr<FaceEle> initial_mass_ele(
new FaceEle(
m_field));
592 post_proc->generateReferenceElementMesh();
603 int main(
int argc,
char *argv[]) {
607 moab::Core mb_instance;
608 moab::Interface &moab = mb_instance;
610 DMType dm_name =
"DMMOFEM";
MoFEM::FaceElementForcesAndSourcesCore FaceEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
#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 DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
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
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
const double v
phase velocity of light in medium (cm/ns)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
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)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
OpCalculateScalarFieldValuesDot OpCalculateScalarValuesDot
auto createTS(MPI_Comm comm)
OpSetContravariantPiolaTransformOnEdge2D OpSetContrariantPiolaTransformOnEdge
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
constexpr double t
plate stiffness
double operator()(const double u, const double v) const
double operator()(const double u, const double v) const
double operator()(const double x, const double y, const double z, const double t) const
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.
Deprecated interface functions.
Section manager is used to create indexes and sections.
Interface for managing meshsets containing materials and boundary conditions.
Get vector field for H-div approximation.
Calculate divergence of vector field.
Get value at integration points for scalar field.
Make Hdiv space from Hcurl space in 2d.
keeps basic data about problem
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
MoFEMErrorCode addDataField(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 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 addBoundaryField(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 boundary.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getDomainFEName() const
Get the Domain FE Name.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode update_slow_rhs(std::string mass_fiedl, boost::shared_ptr< VectorDouble > &mass_ptr)
boost::shared_ptr< MatrixDouble > flux_values_ptr2
boost::shared_ptr< FaceEle > vol_ele_stiff_lhs
std::map< int, BlockData > material_blocks
MoFEMErrorCode update_stiff_rhs()
MoFEMErrorCode apply_BC(std::string flux_field)
boost::shared_ptr< Monitor > monitor_ptr
MoFEMErrorCode push_slow_rhs()
boost::shared_ptr< PreviousData > data1
MoFEMErrorCode update_vol_fe(boost::shared_ptr< FaceEle > &vol_ele, boost::shared_ptr< PreviousData > &data)
Range essential_bdry_ents
MoFEMErrorCode post_proc_fields()
boost::shared_ptr< ForcesAndSourcesCore > null
boost::shared_ptr< VectorDouble > flux_divs_ptr2
MoFEM::Interface & m_field
Simple * simple_interface
MoFEMErrorCode run_analysis()
boost::shared_ptr< VectorDouble > flux_divs_ptr1
MoFEMErrorCode extract_bd_ents(std::string ESSENTIAL, std::string NATURAL)
boost::shared_ptr< BoundaryEle > natural_bdry_ele_slow_rhs
MoFEMErrorCode set_blockData(std::map< int, BlockData > &block_data_map)
RDProblem(MoFEM::Core &core, const int order)
MoFEMErrorCode push_stiff_lhs()
MoFEMErrorCode push_stiff_rhs()
MoFEMErrorCode extract_initial_ents(int block_id, Range &surface)
boost::shared_ptr< VectorDouble > mass_values_ptr2
MoFEMErrorCode apply_IC(std::string mass_field, Range &surface, boost::shared_ptr< FaceEle > &initial_ele, double &init_val)
boost::shared_ptr< VectorDouble > mass_dots_ptr1
boost::shared_ptr< VectorDouble > mass_values_ptr1
boost::shared_ptr< VectorDouble > mass_dots_ptr2
boost::shared_ptr< PreviousData > data2
MoFEMErrorCode setup_system()
MoFEMErrorCode output_result()
std::vector< boost::shared_ptr< PreviousData > > data
boost::shared_ptr< FaceEle > vol_ele_slow_rhs
MoFEMErrorCode update_stiff_lhs()
MoFEMErrorCode set_integration_rule()
boost::shared_ptr< PostProcFaceOnRefinedMesh > post_proc
boost::shared_ptr< MatrixDouble > flux_values_ptr1
boost::shared_ptr< FaceEle > vol_ele_stiff_rhs
double operator()(const double u, const double v, const double dt) const
double operator()(const double u, const double v) const
double operator()(const double u, const double v) const