libMesh
subdomains_ex3.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // <h1>Subdomains Example 3 - Integrating discontinuous data that cuts the mesh</h1>
21 // \author Benjamin S. Kirk
22 // \date 2013
23 
24 
25 // C++ include files that we need
26 #include <iostream>
27 #include <algorithm>
28 #include <math.h>
29 
30 // Basic include file needed for the mesh functionality.
31 #include "libmesh/libmesh.h"
32 #include "libmesh/mesh.h"
33 #include "libmesh/mesh_generation.h"
34 #include "libmesh/mesh_refinement.h"
35 #include "libmesh/quadrature_gauss.h"
36 #include "libmesh/quadrature_composite.h"
37 #include "libmesh/fe.h"
38 #include "libmesh/elem.h"
39 #include "libmesh/parallel.h"
40 #include "libmesh/getpot.h"
41 
42 // Bring in everything from the libMesh namespace
43 using namespace libMesh;
44 
45 // declare the functions we will use
46 void integrate_function (const MeshBase & mesh);
47 
48 // signed distance function
49 const Real radius = 0.5;
50 
51 Real distance (const Point & p)
52 {
53  Point cent(0.8, 0.9);
54  return ((p-cent).norm() - radius);
55 }
56 
57 Real integrand (const Point & p)
58 {
59  return (distance(p) < 0) ? 10. : 1.;
60 }
61 
62 
63 
64 // Begin the main program.
65 int main (int argc, char ** argv)
66 {
67  // Initialize libMesh and any dependent libraries, like in example 2.
68  LibMeshInit init (argc, argv);
69 
70  // This example requires Adaptive Mesh Refinement support - although
71  // it only refines uniformly, the refinement code used is the same
72  // underneath. It also requires libmesh support for Triangle and
73  // Tetgen, which means that libmesh must be configured with
74  // --disable-strict-lgpl for this example to run.
75 #if !defined(LIBMESH_HAVE_TRIANGLE) || !defined(LIBMESH_HAVE_TETGEN) || !defined(LIBMESH_ENABLE_AMR)
76  libmesh_example_requires(false, "--disable-strict-lgpl --enable-amr");
77 #else
78 
79  // Skip this 2D example if libMesh was compiled as 1D-only.
80  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
81 
82  // Read the mesh from file. This is the coarse mesh that will be used
83  // in example 10 to demonstrate adaptive mesh refinement. Here we will
84  // simply read it in and uniformly refine it 5 times before we compute
85  // with it.
86  Mesh mesh(init.comm());
87 
88  {
89  const unsigned int dim =
91 
92  if (dim == 3)
93  std::cout << "Running in 3D" << std::endl;
94 
95  mesh.read ((dim==2) ? "mesh.xda" : "hybrid_3d.xda");
96  }
97 
98  // Create a MeshRefinement object to handle refinement of our mesh.
99  // This class handles all the details of mesh refinement and coarsening.
100  MeshRefinement mesh_refinement (mesh);
101 
102  // Uniformly refine the mesh, by default 4 times.
103  const int n_refinements =
104  libMesh::command_line_next("-n_refinements", 4);
105 
106  mesh_refinement.uniformly_refine (n_refinements);
107 
108  // Print information about the mesh to the screen.
109  mesh.print_info();
110 
111  // integrate the desired function
113 
114  // All done.
115  return 0;
116 #endif
117 }
118 
119 
120 
122 {
123 #if defined(LIBMESH_HAVE_TRIANGLE) && defined(LIBMESH_HAVE_TETGEN)
124 
125  std::vector<Real> vertex_distance;
126 
128  //QGauss qrule (mesh.mesh_dimension(), FIRST);
129 
130  std::unique_ptr<FEBase> fe (FEBase::build (mesh.mesh_dimension(), FEType (FIRST, LAGRANGE)));
131 
132  Real int_val=0.;
133 
134  const std::vector<Point> & q_points = fe->get_xyz();
135  const std::vector<Real> & JxW = fe->get_JxW();
136 
137  for (const auto & elem : mesh.active_local_element_ptr_range())
138  {
139  vertex_distance.clear();
140 
141  for (unsigned int v=0; v<elem->n_vertices(); v++)
142  vertex_distance.push_back (distance(elem->point(v)));
143 
144  qrule.init (*elem, vertex_distance);
145 
146  fe->reinit (elem,
147  &(qrule.get_points()),
148  &(qrule.get_weights()));
149 
150 
151  // TODO: would it be valuable to have the composite quadrature rule sort
152  // from smallest to largest JxW value to help prevent
153  // ... large + small + large + large + small ...
154  // type truncation errors?
155  for (std::size_t qp=0; qp<q_points.size(); qp++)
156  int_val += JxW[qp] * integrand(q_points[qp]);
157  }
158 
159  mesh.comm().sum (int_val);
160 
161  libMesh::out << "\n***********************************\n"
162  << " int_val = " << int_val << std::endl
163  << " exact_val = " << 1*(2*2 - radius*radius*pi) + 10.*(radius*radius*pi)
164  << "\n***********************************\n"
165  << std::endl;
166 #else
168 #endif
169 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
T command_line_next(std::string name, T default_value)
Use GetPot&#39;s search()/next() functions to get following arguments from the command line...
Definition: libmesh.C:1078
virtual void read(const std::string &name, void *mesh_data=nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
const Real radius
unsigned int dim
void sum(T &r) const
MeshBase & mesh
const Parallel::Communicator & comm() const
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:90
virtual void init(const Elem &elem, const std::vector< Real > &vertex_distance_func, unsigned int p_level=0) override
Overrides the base class init() function, and uses the ElemCutter to subdivide the element into "insi...
The libMesh namespace provides an interface to certain functionality in the library.
Real distance(const Point &p)
Real integrand(const Point &p)
This is the MeshBase class.
Definition: mesh_base.h:75
void libmesh_ignore(const Args &...)
void integrate_function(const MeshBase &mesh)
Implements (adaptive) mesh refinement algorithms for a MeshBase.
void print_info(std::ostream &os=libMesh::out, const unsigned int verbosity=0, const bool global=true) const
Prints relevant information about the mesh.
Definition: mesh_base.C:1562
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
auto norm(const T &a) -> decltype(std::abs(a))
Definition: tensor_tools.h:74
This class implements generic composite quadrature rules.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
int main(int argc, char **argv)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
const Real pi
.
Definition: libmesh.h:299
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.