82using namespace common;
84int main(
int argc,
char *argv[])
87 cout <<
"\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this miniapp\n"
88 <<
"is NOT supported with the GPU version of hypre.\n\n";
89 return MFEM_SKIP_RETURN_VALUE;
98 const char *mesh_file =
"../../data/inline-quad.mesh";
100 bool visualization =
true;
101 int ser_ref_levels = 0;
102 int dirichlet_level_set_type = -1;
103 int neumann_level_set_type = -1;
104 bool dirichlet_combo =
false;
108 bool include_cut_cell =
false;
111 args.
AddOption(&mesh_file,
"-m",
"--mesh",
112 "Mesh file to use.");
114 "Finite element order (polynomial degree) or -1 for"
115 " isoparametric space.");
116 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
117 "--no-visualization",
118 "Enable or disable GLVis visualization.");
119 args.
AddOption(&ser_ref_levels,
"-rs",
"--refine-serial",
120 "Number of times to refine the mesh uniformly in serial.");
121 args.
AddOption(&dirichlet_level_set_type,
"-lst",
"--level-set-type",
123 args.
AddOption(&neumann_level_set_type,
"-nlst",
"--neumann-level-set-type",
124 "neumann-level-set-type.");
125 args.
AddOption(&dirichlet_combo,
"-dc",
"--dcombo",
126 "no-dc",
"--no-dcombo",
127 "Combination of two Dirichlet level sets.");
128 args.
AddOption(&ho_terms,
"-ho",
"--high-order",
129 "Additional high-order terms to include");
131 "Nitsche penalty parameter (~1 for 2D, ~10 for 3D).");
132 args.
AddOption(&include_cut_cell,
"-cut",
"--cut",
"-no-cut-cell",
134 "Include or not include elements cut by true boundary.");
135 args.
AddOption(&visport,
"-p",
"--send-port",
"Socket for GLVis.");
148 if (dirichlet_level_set_type <= 0 && neumann_level_set_type <= 0)
150 dirichlet_level_set_type = 1;
152 MFEM_VERIFY((neumann_level_set_type > 0 && ho_terms < 1) ==
false,
153 "Shifted Neumann BC requires extra terms, i.e., -ho >= 1.");
158 if (myid == 0) { device.
Print(); }
161 Mesh mesh(mesh_file, 1, 1);
166 std::cout <<
"Number of elements: " << mesh.
GetNE() << std::endl;
168 MFEM_VERIFY(mesh.
Conforming(),
"AMR capability is not implemented yet!");
171 ParMesh pmesh(MPI_COMM_WORLD, mesh);
176 if (order < 1) { order = 1; }
190 int nodes_cnt = vxyz.
Size()/
dim;
191 if (dirichlet_level_set_type == 3)
193 for (
int i = 0; i < nodes_cnt; i++)
196 vxyz(i+nodes_cnt) = (1.+1.e-4)*vxyz(i+nodes_cnt)-1.e-4;
204 cout <<
"Number of finite element unknowns: " << gtsize << endl;
224 if (dirichlet_level_set_type > 0)
229 filter.
Filter(*dirichlet_dist_coef, level_set_gf);
243 MFEM_VERIFY(dirichlet_level_set_type == 5,
244 "The combo level set example has been only set for"
245 " dirichlet_level_set_type == 5.");
254 if (neumann_level_set_type > 0)
269 for (
int i = 0; i < elem_marker_gf.
Size(); i++)
271 elem_marker_gf(i) = (
real_t)elem_marker[i];
277 "Element Flags", 0, 0, s, s,
"Rjmpc");
289 for (
int i = 0; i < sb_dofs.
Size(); i++)
291 face_dofs(sb_dofs[i]) = 1.0;
297 "Shifted Face Dofs", 0, s, s, s,
"Rjmplo");
315 if (dirichlet_level_set_type == 1 || dirichlet_level_set_type == 2 ||
316 dirichlet_level_set_type == 3)
328 filter.
Filter(combo_dist_coef, filt_gf);
337 "Input Level Set", 0, 2*s, s, s,
"Rjmm");
354 "Distance Vector", s, s, s, s,
"Rjmmpcvv", 1);
362 bool inactive_elements =
false;
363 for (
int i = 0; i < pmesh.
GetNE(); i++)
365 if (!include_cut_cell &&
366 (elem_marker[i] == ShiftedFaceMarker::SBElementType::OUTSIDE ||
367 elem_marker[i] >= ShiftedFaceMarker::SBElementType::CUT))
370 inactive_elements =
true;
372 if (include_cut_cell &&
373 elem_marker[i] == ShiftedFaceMarker::SBElementType::OUTSIDE)
376 inactive_elements =
true;
379 bool inactive_elements_global;
380 MPI_Allreduce(&inactive_elements, &inactive_elements_global, 1,
381 MFEM_MPI_CXX_BOOL, MPI_LOR, MPI_COMM_WORLD);
382 if (inactive_elements_global) { ess_elem.
Append(0); }
389 if (dirichlet_level_set_type == 1 || dirichlet_level_set_type == 4 ||
390 dirichlet_level_set_type == 5 || dirichlet_level_set_type == 6 ||
391 dirichlet_level_set_type == 8 ||
392 neumann_level_set_type == 1 || neumann_level_set_type == 7)
396 else if (dirichlet_level_set_type == 2 || neumann_level_set_type == 2)
400 else if (dirichlet_level_set_type == 3)
404 else { MFEM_ABORT(
"RHS function not set for level set type.\n"); }
411 if (dirichlet_level_set_type == 1 || dirichlet_level_set_type >= 4)
416 else if (dirichlet_level_set_type == 2)
421 else if (dirichlet_level_set_type == 3)
436 if (neumann_level_set_type == 1)
441 else if (neumann_level_set_type == 2)
447 else if (neumann_level_set_type == 7)
452 else if (neumann_level_set_type > 0)
454 MFEM_ABORT(
" Normal vector coefficient not implemented for level set.");
461 int ls_cut_marker = ShiftedFaceMarker::SBElementType::CUT;
464 Array<int> bf_dirichlet_marker(0), bf_neumann_marker(0);
466 if (dirichlet_level_set_type > 0)
481 bf_dirichlet_marker.Append(ls_cut_marker);
493 bf_dirichlet_marker.Append(ls_cut_marker);
499 if (neumann_level_set_type > 0)
501 MFEM_VERIFY(!include_cut_cell,
"include_cut_cell option must be set to"
502 " false for Neumann boundary conditions.");
504 &pmesh, *nbcCoef, *dist_vec,
505 *normalbcCoef, elem_marker,
506 ho_terms, include_cut_cell,
508 bf_neumann_marker.
Append(ls_cut_marker);
520 if (dirichlet_level_set_type > 0)
532 ho_terms), ess_shift_bdr);
536 if (neumann_level_set_type > 0)
561 a.FormLinearSystem(ess_tdof_list, x,
b, A, X, B);
573 a.RecoverFEMSolution(X,
b, x);
576 ofstream mesh_ofs(
"diffusion.mesh");
577 mesh_ofs.precision(8);
579 ofstream sol_ofs(
"diffusion.gf");
580 sol_ofs.precision(8);
600 "Solution", s, 0, s, s,
"Rj");
604 if (dirichlet_level_set_type == 2 || dirichlet_level_set_type == 3 ||
605 (dirichlet_level_set_type == -1 && neumann_level_set_type == 2))
610 for (
int i = 0; i < nodes_cnt; i++)
613 pxyz(1) = vxyz(i+nodes_cnt);
615 if (dirichlet_level_set_type == 2 || neumann_level_set_type == 2)
619 else if (dirichlet_level_set_type == 3)
623 error(i) = std::fabs(x(i) - exact_val);
632 "Error", 2*s, 0, s, s,
"Rj");
638 std::cout <<
"Global L2 error: " << global_error << endl;
643 if (myid == 0) { std::cout << setprecision(10) << norm << std::endl; }
654 delete neumann_dist_coef;
655 delete dirichlet_dist_coef;
656 delete dirichlet_dist_coef_2;
Combination of level sets: +1 inside the true domain, -1 outside.
void Add_Level_Set_Coefficient(Dist_Level_Set_Coefficient &dls_)
Level set coefficient: +1 inside the true domain, -1 outside.
Distance vector to the zero level-set.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
int Size() const
Return the logical size of the array.
int Append(const T &el)
Append element 'el' to array, resize if necessary.
void Mult(const Vector &b, Vector &x) const override
Iterative solution of the linear system using the BiCGSTAB method.
void SetOperator(const Operator &op) override
Set/update the solver for the given operator.
A coefficient that is constant across space and time.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetTime(real_t t)
Set physical time (for time-dependent simulations)
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Class for domain integration .
A general function coefficient.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Arbitrary order H1-conforming (continuous) finite elements.
The BoomerAMG solver in hypre.
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
void SetRelTol(real_t rtol)
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
void SetMaxIter(int max_it)
Arbitrary order "L2-conforming" discontinuous finite elements.
bool Conforming() const
Return a bool indicating whether this mesh is conforming.
void SetAttribute(int i, int attr)
Set the attribute of element i.
void Clear()
Clear the contents of the Mesh.
int GetNE() const
Returns number of elements.
int Dimension() const
Dimension of the reference space used within the elements.
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
void GetNodes(Vector &node_coord) const
void SetNodes(const Vector &node_coord)
Updates the vertex/node locations. Invokes NodesUpdated().
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Array< int > attributes
A list of all unique element attributes used by the Mesh.
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
Pointer to an Operator of a specified type.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
bool Good() const
Return true if the command line options were parsed successfully.
Abstract parallel finite element space.
void ExchangeFaceNbrData()
HYPRE_BigInt GlobalTrueVSize() const
Class for parallel grid function.
void ExchangeFaceNbrData()
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
Returns ||u_ex - u_h||_L2 in parallel for H1 or L2 elements.
void ProjectDiscCoefficient(VectorCoefficient &coeff) override
Project a discontinuous vector coefficient as a grid function on a continuous finite element space....
MFEM_DEPRECATED real_t ComputeL1Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const override
Returns ||u_ex - u_h||_L1 in parallel for H1 or L2 elements.
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void SaveAsOne(const char *fname, int precision=16) const
Class for parallel meshes.
void SetAttributes() override
Determine the sets of unique attribute values in domain and boundary elements.
void SetNodalFESpace(FiniteElementSpace *nfes) override
void PrintAsOne(std::ostream &out=mfem::out, const std::string &comments="") const
Helper class for ParaView visualization data.
void SetLevelsOfDetail(int levels_of_detail_)
void ListShiftedFaceDofs(const Array< int > &elem_marker, Array< int > &sface_dof_list) const
void MarkElements(const ParGridFunction &ls_func, Array< int > &elem_marker)
void ListEssentialTDofs(const Array< int > &elem_marker, const Array< int > &sface_dof_list, Array< int > &ess_tdof_list, Array< int > &ess_shift_bdr) const
Base class for vector Coefficients that optionally depend on time and space.
Vector coefficient defined by a vector GridFunction.
int Size() const
Returns the size of the vector.
IterativeSolver::PrintLevel print_level
virtual void ComputeVectorDistance(Coefficient &zero_level_set, ParGridFunction &distance)
void Filter(ParGridFunction &func, ParGridFunction &ffield)
real_t AvgElementSize(ParMesh &pmesh)
void VisualizeField(socketstream &sock, const char *vishost, int visport, GridFunction &gf, const char *title, int x, int y, int w, int h, const char *keys, bool vec)
real_t rhs_fun_xy_exponent(const Vector &x)
real_t dirichlet_velocity_xy_exponent(const Vector &x)
real_t dirichlet_velocity_xy_sinusoidal(const Vector &x)
real_t rhs_fun_xy_sinusoidal(const Vector &x)
real_t homogeneous(const Vector &x)
Boundary conditions - Dirichlet.
real_t traction_xy_exponent(const Vector &x)
Neumann condition for exponent based solution.
void normal_vector_1(const Vector &x, Vector &p)
real_t rhs_fun_circle(const Vector &x)
f for the Poisson problem (-nabla^2 u = f).
void normal_vector_2(const Vector &x, Vector &p)
Normal vector for level_set_type = 7. Circle centered at [0.5 , 0.6].
PrintLevel & FirstAndLast()