libMesh
mesh_function_dfem.C
Go to the documentation of this file.
1 // Ignore unused parameter warnings coming from cppunit headers
2 #include <libmesh/ignore_warnings.h>
3 #include <cppunit/extensions/HelperMacros.h>
4 #include <cppunit/TestCase.h>
5 #include <libmesh/restore_warnings.h>
6 
7 #include <libmesh/equation_systems.h>
8 #include <libmesh/replicated_mesh.h>
9 #include <libmesh/mesh_generation.h>
10 #include <libmesh/edge_edge2.h>
11 #include <libmesh/face_quad4.h>
12 #include <libmesh/face_tri3.h>
13 #include <libmesh/cell_hex8.h>
14 #include <libmesh/dof_map.h>
15 #include <libmesh/linear_implicit_system.h>
16 #include <libmesh/mesh_refinement.h>
17 #include <libmesh/mesh_function.h>
18 #include <libmesh/numeric_vector.h>
19 
20 #include "test_comm.h"
21 
22 // THE CPPUNIT_TEST_SUITE_END macro expands to code that involves
23 // std::auto_ptr, which in turn produces -Wdeprecated-declarations
24 // warnings. These can be ignored in GCC as long as we wrap the
25 // offending code in appropriate pragmas. We can't get away with a
26 // single ignore_warnings.h inclusion at the beginning of this file,
27 // since the libmesh headers pull in a restore_warnings.h at some
28 // point. We also don't bother restoring warnings at the end of this
29 // file since it's not a header.
30 #include <libmesh/ignore_warnings.h>
31 
32 using namespace libMesh;
33 
35  const Parameters&,
36  const std::string&,
37  const std::string&)
38 {
39  if (p(1) > 0.0)
40  return 0.0;
41  else
42  return 1.0;
43 }
44 
46  const Parameters&,
47  const std::string&,
48  const std::string&)
49 {
50  if (p(1) > 0.0)
51  return 2.0 * p(1) + 1.0;
52  else
53  return p(1) + 1.0;
54 }
55 
56 class MeshfunctionDFEM : public CppUnit::TestCase
57 {
62 public:
63  CPPUNIT_TEST_SUITE( MeshfunctionDFEM );
64 
65  CPPUNIT_TEST( test_point_locator_dfem );
66 
67  CPPUNIT_TEST( test_mesh_function_dfem );
68 
69  CPPUNIT_TEST( test_mesh_function_dfem_grad );
70 
71  CPPUNIT_TEST_SUITE_END();
72 
73 protected:
76 
77  void build_mesh()
78  {
79  _mesh = new ReplicatedMesh(*TestCommWorld);
80 
81  // (0,1) (1,1)
82  // x---------------x
83  // | |
84  // | |
85  // | |
86  // | |
87  // | |
88  // x---------------x
89  // (0,0) (1,0)
90  // | |
91  // | |
92  // | |
93  // | |
94  // x---------------x
95  // (0,-1) (1,-1)
96 
97  _mesh->set_mesh_dimension(2);
98 
99  _mesh->add_point( Point(0.0,-1.0), 4 );
100  _mesh->add_point( Point(1.0,-1.0), 5 );
101  _mesh->add_point( Point(1.0, 0.0), 1 );
102  _mesh->add_point( Point(1.0, 1.0), 2 );
103  _mesh->add_point( Point(0.0, 1.0), 3 );
104  _mesh->add_point( Point(0.0, 0.0), 0 );
105 
106  {
107  Elem* elem_top = _mesh->add_elem( new Quad4 );
108  elem_top->set_node(0) = _mesh->node_ptr(0);
109  elem_top->set_node(1) = _mesh->node_ptr(1);
110  elem_top->set_node(2) = _mesh->node_ptr(2);
111  elem_top->set_node(3) = _mesh->node_ptr(3);
112 
113  Elem* elem_bottom = _mesh->add_elem( new Quad4 );
114  elem_bottom->set_node(0) = _mesh->node_ptr(4);
115  elem_bottom->set_node(1) = _mesh->node_ptr(5);
116  elem_bottom->set_node(2) = _mesh->node_ptr(1);
117  elem_bottom->set_node(3) = _mesh->node_ptr(0);
118  }
119 
120  // libMesh will renumber, but we numbered according to its scheme
121  // anyway. We do this because when we call uniformly_refine subsequently,
122  // it's going use skip_renumber=false.
123  _mesh->prepare_for_use(false /*skip_renumber*/);
124 
125  // get a point locator
126  _point_locator = _mesh->sub_point_locator();
127  }
128 
129 public:
130  void setUp()
131  {
132  this->build_mesh();
133  }
134 
135  void tearDown()
136  {
137  delete _mesh;
138  }
139 
140  // test that point locator works correctly
142  {
143  // this point is in bottom element only
144  Point interior(0.5, -0.5, 0.0);
145 
146  // this point is on the face between bottom and top
147  Point face(0.5, 0.0, 0.0);
148 
149  // test interior point
150  std::set<const Elem *> int_cand;
151  (*_point_locator)(interior, int_cand, libmesh_nullptr);
152  CPPUNIT_ASSERT (int_cand.size() == 1);
153  for (std::set<const Elem *>::iterator it = int_cand.begin(); it != int_cand.end(); ++it)
154  CPPUNIT_ASSERT ((*it)->id() == 1);
155 
156  // test interior point
157  std::set<const Elem *> face_cand;
158  (*_point_locator)(face, face_cand, libmesh_nullptr);
159  CPPUNIT_ASSERT (face_cand.size() == 2);
160  int array[2] = {0, 0};
161  for (std::set<const Elem *>::iterator it = face_cand.begin(); it != face_cand.end(); ++it)
162  {
163  // first test that array entry hasn't been set before
164  CPPUNIT_ASSERT (array[(*it)->id()] == 0);
165  array[(*it)->id()] = 1;
166  }
167  CPPUNIT_ASSERT (array[0] == 1);
168  CPPUNIT_ASSERT (array[1] == 1);
169  }
170 
171  // test that mesh function works correctly
173  {
174  // this point is in top element only
175  Point top(0.5, 0.5, 0.0);
176 
177  // this point is in bottom element only
178  Point bottom(0.5, -0.5, 0.0);
179 
180  // this point is on the face between bottom and top
181  Point face(0.5, 0.0, 0.0);
182 
183  // set up some necessary objects
184  EquationSystems es(*_mesh);
185  System & sys = es.add_system<System> ("SimpleSystem");
186  sys.add_variable("u", CONSTANT, MONOMIAL);
187 
188  es.init();
190 
191  std::vector<unsigned int> variables;
192  sys.get_all_variable_numbers(variables);
193  UniquePtr<NumericVector<Number>> mesh_function_vector =
195  mesh_function_vector->init(sys.n_dofs(), false, SERIAL);
196  sys.solution->localize(*mesh_function_vector);
197 
198  MeshFunction mesh_function (es, *mesh_function_vector,
199  sys.get_dof_map(), variables);
200  mesh_function.init();
201 
202  // test mesh function in top
203  std::map<const Elem *, Number> top_val = mesh_function.discontinuous_value(top);
204 
205  // check that there is only one value
206  CPPUNIT_ASSERT (top_val.size() == 1);
207 
208  // check that this one value is the right one
209  for (std::map<const Elem *, Number>::const_iterator it = top_val.begin(); it != top_val.end(); ++it)
210  CPPUNIT_ASSERT (it->first->id() == 0 && std::abs(it->second) < 1.0e-10);
211 
212  // test mesh function in bottom
213  std::map<const Elem *, Number> bottom_val = mesh_function.discontinuous_value(bottom);
214 
215  // check that there is only one value
216  CPPUNIT_ASSERT (bottom_val.size() == 1);
217 
218  // check that this one value is the right one
219  for (std::map<const Elem *, Number>::const_iterator it = bottom_val.begin(); it != bottom_val.end(); ++it)
220  CPPUNIT_ASSERT (it->first->id() == 1 && std::abs(it->second - 1.0) < 1.0e-10);
221 
222  // test mesh function in face
223  std::map<const Elem *, Number> face_val = mesh_function.discontinuous_value(face);
224 
225  // check that there are two values
226  CPPUNIT_ASSERT (face_val.size() == 2);
227 
228  // check that the two values are attached to the right element
229  for (std::map<const Elem *, Number>::const_iterator it = face_val.begin(); it != face_val.end(); ++it)
230  if (it->first->id() == 0)
231  CPPUNIT_ASSERT (std::abs(it->second) < 1.0e-10);
232  else if (it->first->id() == 1)
233  CPPUNIT_ASSERT (std::abs(it->second - 1.0) < 1.0e-10);
234  else
235  CPPUNIT_ASSERT (false);
236 
237  }
238 
239  // test that gradient function works correctly
241  {
242  // this point is in top element only
243  Point top(0.5, 0.5, 0.0);
244 
245  // this point is in bottom element only
246  Point bottom(0.5, -0.5, 0.0);
247 
248  // this point is on the face between bottom and top
249  Point face(0.5, 0.0, 0.0);
250 
251  // set up some necessary objects
252  EquationSystems es(*_mesh);
253  System & sys = es.add_system<System> ("SimpleSystem");
254  sys.add_variable("u", FIRST, LAGRANGE);
255 
256  es.init();
258 
259  std::vector<unsigned int> variables;
260  sys.get_all_variable_numbers(variables);
261  UniquePtr<NumericVector<Number>> mesh_function_vector =
263  mesh_function_vector->init(sys.n_dofs(), false, SERIAL);
264  sys.solution->localize(*mesh_function_vector);
265 
266  MeshFunction mesh_function (es, *mesh_function_vector,
267  sys.get_dof_map(), variables);
268  mesh_function.init();
269 
270  // test mesh function in top
271  std::map<const Elem *, Gradient> top_val = mesh_function.discontinuous_gradient(top);
272 
273  // check that there is only one value
274  CPPUNIT_ASSERT (top_val.size() == 1);
275 
276  // check that this one value is the right one
277  for (std::map<const Elem *, Gradient>::const_iterator it = top_val.begin(); it != top_val.end(); ++it)
278  CPPUNIT_ASSERT (it->first->id() == 0 && std::abs(it->second(1) - 2.0) < 1.0e-10);
279 
280  // test mesh function in bottom
281  std::map<const Elem *, Gradient> bottom_val = mesh_function.discontinuous_gradient(bottom);
282 
283  // check that there is only one value
284  CPPUNIT_ASSERT (bottom_val.size() == 1);
285 
286  // check that this one value is the right one
287  for (std::map<const Elem *, Gradient>::const_iterator it = bottom_val.begin(); it != bottom_val.end(); ++it)
288  CPPUNIT_ASSERT (it->first->id() == 1 && std::abs(it->second(1) - 1.0) < 1.0e-10);
289 
290  // test mesh function in face
291  std::map<const Elem *, Gradient> face_val = mesh_function.discontinuous_gradient(face);
292 
293  // check that there are two values
294  CPPUNIT_ASSERT (face_val.size() == 2);
295 
296  // check that the two values are attached to the right element
297  for (std::map<const Elem *, Gradient>::const_iterator it = face_val.begin(); it != face_val.end(); ++it)
298  if (it->first->id() == 0)
299  CPPUNIT_ASSERT (std::abs(it->second(1) - 2.0) < 1.0e-10);
300  else if (it->first->id() == 1)
301  CPPUNIT_ASSERT (std::abs(it->second(1) - 1.0) < 1.0e-10);
302  else
303  CPPUNIT_ASSERT (false);
304  }
305 };
306 
307 
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id) libmesh_override
functions for adding /deleting nodes elements.
Number position_function2(const Point &p, const Parameters &, const std::string &, const std::string &)
double abs(double a)
This is the EquationSystems class.
virtual const Node * node_ptr(const dof_id_type i) const libmesh_override
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1941
This class provides the ability to map between arbitrary, user-defined strings and several data types...
Definition: parameters.h:63
The QUAD4 is an element in 2D composed of 4 nodes.
Definition: face_quad4.h:49
virtual Elem * add_elem(Elem *e) libmesh_override
Add elem e to the end of the element array.
UniquePtr< PointLocatorBase > _point_locator
libMesh::Parallel::Communicator * TestCommWorld
Definition: driver.C:28
ImplicitSystem & sys
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
const class libmesh_nullptr_t libmesh_nullptr
void project_solution(FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr) const
Projects arbitrary functions onto the current solution.
The libMesh namespace provides an interface to certain functionality in the library.
void test_mesh_function_dfem_grad()
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
const DofMap & get_dof_map() const
Definition: system.h:2030
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
void prepare_for_use(const bool skip_renumber_nodes_and_elements=false, const bool skip_find_neighbors=false)
Prepare a newly created (or read) mesh for use.
Definition: mesh_base.C:174
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
void set_mesh_dimension(unsigned char d)
Resets the logical dimension of the mesh.
Definition: mesh_base.h:199
Number position_function(const Point &p, const Parameters &, const std::string &, const std::string &)
const Parallel::Communicator & comm() const
Parameters parameters
Data structure holding arbitrary parameters.
CPPUNIT_TEST_SUITE_REGISTRATION(MeshfunctionDFEM)
dof_id_type n_dofs() const
Definition: system.C:148
ReplicatedMesh * _mesh
virtual void init()
Initialize all the systems.
This class provides function-like objects for data distributed over a mesh.
Definition: mesh_function.h:53
void get_all_variable_numbers(std::vector< unsigned int > &all_variable_numbers) const
Fills all_variable_numbers with all the variable numbers for the variables that have been added to th...
Definition: system.C:1278
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
UniquePtr< PointLocatorBase > sub_point_locator() const
Definition: mesh_base.C:534
virtual void init() libmesh_override
Override the FunctionBase::init() member function by calling our own and specifying the Trees::NODES ...
Definition: mesh_function.h:96