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 UniquePtr<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 
105  elem_it = _mesh->active_local_elements_begin(),
106  elem_end = _mesh->active_local_elements_end();
107 
108  if (elem_it == elem_end)
109  _elem = libmesh_nullptr;
110  else
111  _elem = *elem_it;
112 
113  _sys->get_dof_map().dof_indices(_elem, _dof_indices);
114 
115  _nx = 10;
116  _ny = _nx * (_dim > 1);
117  _nz = _nx * (_dim > 2);
118  }
119 
120  void tearDown()
121  {
122  delete _fe;
123  delete _es;
124  delete _mesh;
125  }
126 
127  void testU()
128  {
129  // Handle the "more processors than elements" case
130  if (!_elem)
131  return;
132 
133  // These tests require exceptions to be enabled because a
134  // TypeTensor::solve() call down in Elem::contains_point()
135  // actually throws a non-fatal exception for a certain Point which
136  // is not in the Elem. When exceptions are not enabled, this test
137  // simply aborts.
138 #ifdef LIBMESH_ENABLE_EXCEPTIONS
139  for (unsigned int i=0; i != _nx; ++i)
140  for (unsigned int j=0; j != _ny; ++j)
141  for (unsigned int k=0; k != _nz; ++k)
142  {
143  Real x = (Real(i)/_nx), y = 0, z = 0;
144  Point p = x;
145  if (j > 0)
146  p(1) = y = (Real(j)/_ny);
147  if (k > 0)
148  p(2) = z = (Real(k)/_nz);
149  if (!_elem->contains_point(p))
150  continue;
151 
152  std::vector<Point> master_points
153  (1, FEInterface::inverse_map(_dim, _fe->get_fe_type(), _elem, p));
154 
155  _fe->reinit(_elem, &master_points);
156 
157  Number u = 0;
158  for (std::size_t d = 0; d != _dof_indices.size(); ++d)
159  u += _fe->get_phi()[d][0] * (*_sys->current_local_solution)(_dof_indices[d]);
160 
161  CPPUNIT_ASSERT_DOUBLES_EQUAL
162  (libmesh_real(u),
163  libmesh_real(x + 0.25*y + 0.0625*z),
165  }
166 #endif
167  }
168 
169  void testGradU()
170  {
171  // Handle the "more processors than elements" case
172  if (!_elem)
173  return;
174 
175  // These tests require exceptions to be enabled because a
176  // TypeTensor::solve() call down in Elem::contains_point()
177  // actually throws a non-fatal exception for a certain Point which
178  // is not in the Elem. When exceptions are not enabled, this test
179  // simply aborts.
180 #ifdef LIBMESH_ENABLE_EXCEPTIONS
181  for (unsigned int i=0; i != _nx; ++i)
182  for (unsigned int j=0; j != _ny; ++j)
183  for (unsigned int k=0; k != _nz; ++k)
184  {
185  Point p(Real(i)/_nx);
186  if (j > 0)
187  p(1) = Real(j)/_ny;
188  if (k > 0)
189  p(2) = Real(k)/_ny;
190  if (!_elem->contains_point(p))
191  continue;
192 
193  std::vector<Point> master_points
194  (1, FEInterface::inverse_map(_dim, _fe->get_fe_type(), _elem, p));
195 
196  _fe->reinit(_elem, &master_points);
197 
198  Gradient grad_u = 0;
199  for (std::size_t d = 0; d != _dof_indices.size(); ++d)
200  grad_u += _fe->get_dphi()[d][0] * (*_sys->current_local_solution)(_dof_indices[d]);
201 
202  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(grad_u(0)), 1.0,
203  TOLERANCE*sqrt(TOLERANCE));
204  if (_dim > 1)
205  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(grad_u(1)), 0.25,
206  TOLERANCE*sqrt(TOLERANCE));
207  if (_dim > 2)
208  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(grad_u(2)), 0.0625,
209  TOLERANCE*sqrt(TOLERANCE));
210  }
211 #endif
212  }
213 
215  {
216  // Handle the "more processors than elements" case
217  if (!_elem)
218  return;
219 
220  // These tests require exceptions to be enabled because a
221  // TypeTensor::solve() call down in Elem::contains_point()
222  // actually throws a non-fatal exception for a certain Point which
223  // is not in the Elem. When exceptions are not enabled, this test
224  // simply aborts.
225 #ifdef LIBMESH_ENABLE_EXCEPTIONS
226  for (unsigned int i=0; i != _nx; ++i)
227  for (unsigned int j=0; j != _ny; ++j)
228  for (unsigned int k=0; k != _nz; ++k)
229  {
230  Point p(Real(i)/_nx);
231  if (j > 0)
232  p(1) = Real(j)/_ny;
233  if (k > 0)
234  p(2) = Real(k)/_ny;
235  if (!_elem->contains_point(p))
236  continue;
237 
238  std::vector<Point> master_points
239  (1, FEInterface::inverse_map(_dim, _fe->get_fe_type(), _elem, p));
240 
241  _fe->reinit(_elem, &master_points);
242 
243  Number grad_u_x = 0, grad_u_y = 0, grad_u_z = 0;
244  for (std::size_t d = 0; d != _dof_indices.size(); ++d)
245  {
246  grad_u_x += _fe->get_dphidx()[d][0] * (*_sys->current_local_solution)(_dof_indices[d]);
247 #if LIBMESH_DIM > 1
248  grad_u_y += _fe->get_dphidy()[d][0] * (*_sys->current_local_solution)(_dof_indices[d]);
249 #endif
250 #if LIBMESH_DIM > 2
251  grad_u_z += _fe->get_dphidz()[d][0] * (*_sys->current_local_solution)(_dof_indices[d]);
252 #endif
253  }
254 
255  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(grad_u_x), 1.0,
256  TOLERANCE*sqrt(TOLERANCE));
257  if (_dim > 1)
258  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(grad_u_y), 0.25,
259  TOLERANCE*sqrt(TOLERANCE));
260  if (_dim > 2)
261  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(grad_u_z), 0.0625,
262  TOLERANCE*sqrt(TOLERANCE));
263  }
264 #endif
265  }
266 
267 };
268 
269 // THE CPPUNIT_TEST_SUITE_END macro expands to code that involves
270 // std::auto_ptr, which in turn produces -Wdeprecated-declarations
271 // warnings. These can be ignored in GCC as long as we wrap the
272 // offending code in appropriate pragmas. We'll put an
273 // ignore_warnings at the end of this file so it's the last warnings
274 // related header that our including code sees.
275 #include <libmesh/ignore_warnings.h>
276 
277 #define INSTANTIATE_FETEST(order, family, elemtype) \
278  class FETest_##order##_##family##_##elemtype : public FETest<order, family, elemtype> { \
279  public: \
280  CPPUNIT_TEST_SUITE( FETest_##order##_##family##_##elemtype ); \
281  FETEST \
282  CPPUNIT_TEST_SUITE_END(); \
283  }; \
284  \
285  CPPUNIT_TEST_SUITE_REGISTRATION( FETest_##order##_##family##_##elemtype );
286 
287 #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:178
EquationSystems * _es
Definition: fe_test.h:70
virtual element_iterator active_local_elements_begin() libmesh_override
This is the EquationSystems class.
static UniquePtr< Elem > build(const ElemType type, Elem *p=libmesh_nullptr)
Definition: elem.C:238
void tearDown()
Definition: fe_test.h:120
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
const std::vector< std::vector< OutputShape > > & get_dphidz() const
Definition: fe_base.h:256
libMesh::Parallel::Communicator * TestCommWorld
Definition: driver.C:28
The definition of the const_element_iterator struct.
Definition: mesh_base.h:1494
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
void testGradUComp()
Definition: fe_test.h:214
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
UniquePtr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1535
virtual element_iterator active_local_elements_end() libmesh_override
const class libmesh_nullptr_t libmesh_nullptr
static const Real TOLERANCE
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.
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
FEType get_fe_type() const
Definition: fe_abstract.h:421
void setUp()
Definition: fe_test.h:73
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
std::vector< dof_id_type > _dof_indices
Definition: fe_test.h:66
unsigned int _nz
Definition: fe_test.h:64
Number linear_test(const Point &p, const Parameters &, const std::string &, const std::string &)
Definition: fe_test.h:32
PetscErrorCode Vec x
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
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
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.
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)
Definition: fe_interface.C:569
Elem * _elem
Definition: fe_test.h:65
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:216
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const std::vector< std::vector< OutputShape > > & get_dphidx() const
Definition: fe_base.h:240
Parameters parameters
Data structure holding arbitrary parameters.
FEBase * _fe
Definition: fe_test.h:67
virtual void init()
Initialize all the systems.
const std::vector< std::vector< OutputShape > > & get_dphidy() const
Definition: fe_base.h:248
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:50
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
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.
virtual bool contains_point(const Point &p, Real tol=TOLERANCE) const
Definition: elem.C:2448
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr)=0
This is at the core of this class.
void testGradU()
Definition: fe_test.h:169
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.
Definition: dof_map.C:1917
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
void testU()
Definition: fe_test.h:127