libMesh
fe_test.h
Go to the documentation of this file.
1 #ifndef __fe_test_h__
2 #define __fe_test_h__
3 
4 #include "test_comm.h"
5 
6 #include <libmesh/dof_map.h>
7 #include <libmesh/elem.h>
8 #include <libmesh/equation_systems.h>
9 #include <libmesh/fe_base.h>
10 #include <libmesh/fe_interface.h>
11 #include <libmesh/mesh.h>
12 #include <libmesh/mesh_generation.h>
13 #include <libmesh/numeric_vector.h>
14 #include <libmesh/system.h>
15 
16 // Ignore unused parameter warnings coming from cppunit headers
17 #include <libmesh/ignore_warnings.h>
18 #include <cppunit/extensions/HelperMacros.h>
19 #include <cppunit/TestCase.h>
20 #include <libmesh/restore_warnings.h>
21 
22 #include <vector>
23 
24 #define FETEST \
25  CPPUNIT_TEST( testU ); \
26  CPPUNIT_TEST( testGradU ); \
27  CPPUNIT_TEST( testGradUComp );
28 
29 using namespace libMesh;
30 
31 inline
33  const Parameters&,
34  const std::string&,
35  const std::string&)
36 {
37  const Real & x = p(0);
38  const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
39  const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
40 
41  return x + 0.25*y + 0.0625*z;
42 }
43 
44 inline
46  const Parameters&,
47  const std::string&,
48  const std::string&)
49 {
50  Gradient grad = 1;
51  if (LIBMESH_DIM > 1)
52  grad(1) = 0.25;
53  if (LIBMESH_DIM > 2)
54  grad(2) = 0.0625;
55 
56  return grad;
57 }
58 
59 
60 template <Order order, FEFamily family, ElemType elem_type>
61 class FETest : public CppUnit::TestCase {
62 
63 private:
64  unsigned int _dim, _nx, _ny, _nz;
66  std::vector<dof_id_type> _dof_indices;
71 
72 public:
73  void setUp()
74  {
75  _mesh = new Mesh(*TestCommWorld);
76  const std::unique_ptr<Elem> test_elem = Elem::build(elem_type);
77  _dim = test_elem->dim();
78  const unsigned int ny = _dim > 1;
79  const unsigned int nz = _dim > 2;
80 
82  1, ny, nz,
83  0., 1., 0., ny, 0., nz,
84  elem_type);
85 
86  _es = new EquationSystems(*_mesh);
87  _sys = &(_es->add_system<System> ("SimpleSystem"));
88  _sys->add_variable("u", order, family);
89  _es->init();
91 
92  FEType fe_type = _sys->variable_type(0);
93  _fe = FEBase::build(_dim, fe_type).release();
94  _fe->get_phi();
95  _fe->get_dphi();
96  _fe->get_dphidx();
97 #if LIBMESH_DIM > 1
98  _fe->get_dphidy();
99 #endif
100 #if LIBMESH_DIM > 2
101  _fe->get_dphidz();
102 #endif
103 
104  auto rng = _mesh->active_local_element_ptr_range();
105  _elem = rng.begin() == rng.end() ? nullptr : *(rng.begin());
106 
107  _sys->get_dof_map().dof_indices(_elem, _dof_indices);
108 
109  _nx = 10;
110  _ny = _nx * (_dim > 1);
111  _nz = _nx * (_dim > 2);
112  }
113 
114  void tearDown()
115  {
116  delete _fe;
117  delete _es;
118  delete _mesh;
119  }
120 
121  void testU()
122  {
123  // Handle the "more processors than elements" case
124  if (!_elem)
125  return;
126 
127  // These tests require exceptions to be enabled because a
128  // TypeTensor::solve() call down in Elem::contains_point()
129  // actually throws a non-fatal exception for a certain Point which
130  // is not in the Elem. When exceptions are not enabled, this test
131  // simply aborts.
132 #ifdef LIBMESH_ENABLE_EXCEPTIONS
133  for (unsigned int i=0; i != _nx; ++i)
134  for (unsigned int j=0; j != _ny; ++j)
135  for (unsigned int k=0; k != _nz; ++k)
136  {
137  Real x = (Real(i)/_nx), y = 0, z = 0;
138  Point p = x;
139  if (j > 0)
140  p(1) = y = (Real(j)/_ny);
141  if (k > 0)
142  p(2) = z = (Real(k)/_nz);
143  if (!_elem->contains_point(p))
144  continue;
145 
146  std::vector<Point> master_points
147  (1, FEInterface::inverse_map(_dim, _fe->get_fe_type(), _elem, p));
148 
149  _fe->reinit(_elem, &master_points);
150 
151  Number u = 0;
152  for (std::size_t d = 0; d != _dof_indices.size(); ++d)
153  u += _fe->get_phi()[d][0] * (*_sys->current_local_solution)(_dof_indices[d]);
154 
155  CPPUNIT_ASSERT_DOUBLES_EQUAL
156  (libmesh_real(u),
157  libmesh_real(x + 0.25*y + 0.0625*z),
159  }
160 #endif
161  }
162 
163  void testGradU()
164  {
165  // Handle the "more processors than elements" case
166  if (!_elem)
167  return;
168 
169  // These tests require exceptions to be enabled because a
170  // TypeTensor::solve() call down in Elem::contains_point()
171  // actually throws a non-fatal exception for a certain Point which
172  // is not in the Elem. When exceptions are not enabled, this test
173  // simply aborts.
174 #ifdef LIBMESH_ENABLE_EXCEPTIONS
175  for (unsigned int i=0; i != _nx; ++i)
176  for (unsigned int j=0; j != _ny; ++j)
177  for (unsigned int k=0; k != _nz; ++k)
178  {
179  Point p(Real(i)/_nx);
180  if (j > 0)
181  p(1) = Real(j)/_ny;
182  if (k > 0)
183  p(2) = Real(k)/_ny;
184  if (!_elem->contains_point(p))
185  continue;
186 
187  std::vector<Point> master_points
188  (1, FEInterface::inverse_map(_dim, _fe->get_fe_type(), _elem, p));
189 
190  _fe->reinit(_elem, &master_points);
191 
192  Gradient grad_u = 0;
193  for (std::size_t d = 0; d != _dof_indices.size(); ++d)
194  grad_u += _fe->get_dphi()[d][0] * (*_sys->current_local_solution)(_dof_indices[d]);
195 
196  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(grad_u(0)), 1.0,
198  if (_dim > 1)
199  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(grad_u(1)), 0.25,
201  if (_dim > 2)
202  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(grad_u(2)), 0.0625,
204  }
205 #endif
206  }
207 
209  {
210  // Handle the "more processors than elements" case
211  if (!_elem)
212  return;
213 
214  // These tests require exceptions to be enabled because a
215  // TypeTensor::solve() call down in Elem::contains_point()
216  // actually throws a non-fatal exception for a certain Point which
217  // is not in the Elem. When exceptions are not enabled, this test
218  // simply aborts.
219 #ifdef LIBMESH_ENABLE_EXCEPTIONS
220  for (unsigned int i=0; i != _nx; ++i)
221  for (unsigned int j=0; j != _ny; ++j)
222  for (unsigned int k=0; k != _nz; ++k)
223  {
224  Point p(Real(i)/_nx);
225  if (j > 0)
226  p(1) = Real(j)/_ny;
227  if (k > 0)
228  p(2) = Real(k)/_ny;
229  if (!_elem->contains_point(p))
230  continue;
231 
232  std::vector<Point> master_points
233  (1, FEInterface::inverse_map(_dim, _fe->get_fe_type(), _elem, p));
234 
235  _fe->reinit(_elem, &master_points);
236 
237  Number grad_u_x = 0, grad_u_y = 0, grad_u_z = 0;
238  for (std::size_t d = 0; d != _dof_indices.size(); ++d)
239  {
240  grad_u_x += _fe->get_dphidx()[d][0] * (*_sys->current_local_solution)(_dof_indices[d]);
241 #if LIBMESH_DIM > 1
242  grad_u_y += _fe->get_dphidy()[d][0] * (*_sys->current_local_solution)(_dof_indices[d]);
243 #endif
244 #if LIBMESH_DIM > 2
245  grad_u_z += _fe->get_dphidz()[d][0] * (*_sys->current_local_solution)(_dof_indices[d]);
246 #endif
247  }
248 
249  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(grad_u_x), 1.0,
251  if (_dim > 1)
252  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(grad_u_y), 0.25,
254  if (_dim > 2)
255  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(grad_u_z), 0.0625,
257  }
258 #endif
259  }
260 
261 };
262 
263 // THE CPPUNIT_TEST_SUITE_END macro expands to code that involves
264 // std::auto_ptr, which in turn produces -Wdeprecated-declarations
265 // warnings. These can be ignored in GCC as long as we wrap the
266 // offending code in appropriate pragmas. We'll put an
267 // ignore_warnings at the end of this file so it's the last warnings
268 // related header that our including code sees.
269 #include <libmesh/ignore_warnings.h>
270 
271 #define INSTANTIATE_FETEST(order, family, elemtype) \
272  class FETest_##order##_##family##_##elemtype : public FETest<order, family, elemtype> { \
273  public: \
274  CPPUNIT_TEST_SUITE( FETest_##order##_##family##_##elemtype ); \
275  FETEST \
276  CPPUNIT_TEST_SUITE_END(); \
277  }; \
278  \
279  CPPUNIT_TEST_SUITE_REGISTRATION( FETest_##order##_##family##_##elemtype );
280 
281 #endif // #ifdef __fe_test_h__
T libmesh_real(T a)
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:179
EquationSystems * _es
Definition: fe_test.h:70
This is the EquationSystems class.
void tearDown()
Definition: fe_test.h:114
This class provides the ability to map between arbitrary, user-defined strings and several data types...
Definition: parameters.h:63
System * _sys
Definition: fe_test.h:69
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
virtual void init()
Initialize all the systems.
MetaPhysicL::DualNumber< T, D > sqrt(const MetaPhysicL::DualNumber< T, D > &in)
void testGradUComp()
Definition: fe_test.h:208
This is the base class from which all geometric element types are derived.
Definition: elem.h:101
libMesh::Parallel::Communicator * TestCommWorld
static const Real TOLERANCE
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
The libMesh namespace provides an interface to certain functionality in the library.
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
void setUp()
Definition: fe_test.h:73
const std::vector< std::vector< OutputShape > > & get_dphidx() const
Definition: fe_base.h:239
FEType get_fe_type() const
Definition: fe_abstract.h:428
std::vector< dof_id_type > _dof_indices
Definition: fe_test.h:66
unsigned int _nz
Definition: fe_test.h:64
void project_solution(FunctionBase< Number > *f, FunctionBase< Gradient > *g=nullptr) const
Projects arbitrary functions onto the current solution.
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:215
Number linear_test(const Point &p, const Parameters &, const std::string &, const std::string &)
Definition: fe_test.h:32
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:92
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)=0
This is at the core of this class.
const std::vector< std::vector< OutputShape > > & get_dphidy() const
Definition: fe_base.h:247
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Elem * _elem
Definition: fe_test.h:65
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2217
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Parameters parameters
Data structure holding arbitrary parameters.
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1535
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.
virtual SimpleRange< element_iterator > active_local_element_ptr_range() override
FEBase * _fe
Definition: fe_test.h:67
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=nullptr)
Adds the variable var to the list of variables for this system.
Gradient linear_test_grad(const Point &, const Parameters &, const std::string &, const std::string &)
Definition: fe_test.h:45
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:51
const DofMap & get_dof_map() const
Definition: system.h:2083
Mesh * _mesh
Definition: fe_test.h:68
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
Definition: fe_test.h:61
virtual bool contains_point(const Point &p, Real tol=TOLERANCE) const
void build_cube(UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
Builds a (elements) cube.
This class forms the foundation from which generic finite elements may be derived.
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:207
void testGradU()
Definition: fe_test.h:163
const std::vector< std::vector< OutputShape > > & get_dphidz() const
Definition: fe_base.h:255
void testU()
Definition: fe_test.h:121