COOLFluiD  Release kernel
COOLFluiD is a Collaborative Simulation Environment (CSE) focused on complex MultiPhysics simulations.
utest-mesh-lagrangep1-quad2d.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010-2013 von Karman Institute for Fluid Dynamics, Belgium
2 //
3 // This software is distributed under the terms of the
4 // GNU Lesser General Public License version 3 (LGPLv3).
5 // See doc/lgpl.txt and doc/gpl.txt for the license text.
6 
7 #define BOOST_TEST_DYN_LINK
8 #define BOOST_TEST_MODULE "Test module for the ETYPE shapefunction"
9 
10 #include <boost/assign/list_of.hpp>
11 #include <boost/test/unit_test.hpp>
12 
13 #include "common/Log.hpp"
14 
15 #include "common/Table.hpp"
16 #include "mesh/ContinuousDictionary.hpp"
17 #include "mesh/Integrators/Gauss.hpp"
18 #include "mesh/LagrangeP1/Quad2D.hpp"
19 #include "mesh/Elements.hpp"
20 
21 #include "Tools/Testing/Difference.hpp"
22 
23 using namespace boost::assign;
24 using namespace cf3;
25 using namespace cf3::common;
26 using namespace cf3::mesh;
27 using namespace cf3::mesh::Integrators;
28 using namespace cf3::mesh::LagrangeP1;
29 using namespace cf3::Tools::Testing;
30 
32 
33 typedef Quad2D ETYPE;
34 
35 struct LagrangeP1Quad2DFixture
36 {
37  typedef ETYPE::NodesT NodesT;
39  LagrangeP1Quad2DFixture() : mapped_coords(0.1, 0.8), nodes((NodesT() << 0.5, 0.3,
40  1.1, 1.2,
41  1.35, 1.9,
42  0.8, 2.1).finished())
43  {
44  }
45 
47  ~LagrangeP1Quad2DFixture()
48  {
49  }
51 
52  const ETYPE::MappedCoordsT mapped_coords;
53  const NodesT nodes;
54 
55  struct ConstFunctor
56  {
57  ConstFunctor(const NodesT& node_list) : m_nodes(node_list) {}
58 
59  Real operator()() const
60  {
61  return ETYPE::jacobian_determinant(mapped_coords, m_nodes);
62  }
63  ETYPE::MappedCoordsT mapped_coords;
64  private:
65  const NodesT& m_nodes;
66  };
67 };
68 
70 
71 BOOST_FIXTURE_TEST_SUITE( ETYPESuite, LagrangeP1Quad2DFixture )
72 
73 
75 BOOST_AUTO_TEST_CASE( Volume )
76 {
77  ETYPE::NodesT nodes_quad2D;
78  nodes_quad2D <<
79  0.0, 0.0,
80  1.0, 0.0,
81  1.0, 1.0,
82  0.0, 1.0;
83  BOOST_CHECK_EQUAL(ETYPE::volume(nodes_quad2D), 1.0);
84 }
85 
86 BOOST_AUTO_TEST_CASE( Element )
87 {
88  // Create a Elements component
89  boost::shared_ptr<Elements> comp = allocate_component<Elements>("comp");
90  boost::shared_ptr<Dictionary> nodes = allocate_component<ContinuousDictionary>("nodes");
91  comp->initialize("cf3.mesh.LagrangeP1.Quad2D",*nodes);
92  BOOST_CHECK_EQUAL(comp->element_type().shape(), GeoShape::QUAD);
93  BOOST_CHECK_EQUAL(comp->element_type().nb_faces(), (Uint) 4);
94 
95  // Check volume calculation
96  ETYPE::NodesT coord;
97  coord <<
98  15, 15,
99  40, 25,
100  25, 30,
101  30, 40;
102 
103  BOOST_CHECK_EQUAL(ETYPE::volume(coord), 150.);
104 }
105 
106 BOOST_AUTO_TEST_CASE( computeShapeFunction )
107 {
108  const ETYPE::SF::ValueT reference_result(0.045, 0.055, 0.495, 0.405);
109  ETYPE::SF::ValueT result;
110  ETYPE::SF::compute_value(mapped_coords, result);
111  Accumulator accumulator;
112  vector_test(result, reference_result, accumulator);
113  BOOST_CHECK_LT(boost::accumulators::max(accumulator.ulps), 10); // Maximal difference can't be greater than 10 times the least representable unit
114 }
115 
116 BOOST_AUTO_TEST_CASE( computeMappedCoordinates )
117 {
118  Accumulator accumulator;
119  ETYPE::CoordsT test_coords(0.9375, 1.375); // center of the element
120  ETYPE::MappedCoordsT result;
121  ETYPE::compute_mapped_coordinate(test_coords, nodes, result);
122  BOOST_CHECK_LT(std::abs(result[0]), 1e-12);
123  BOOST_CHECK_LT(std::abs(result[1]), 1e-12);// sqrt from the expression gives too many ULPS in difference for Accumulator
124 
125  test_coords = nodes.row(0);
126  ETYPE::compute_mapped_coordinate(test_coords, nodes, result);
127  vector_test(result,ETYPE::MappedCoordsT(-1.0, -1.0),accumulator);
128  BOOST_CHECK_LT(boost::accumulators::max(accumulator.ulps), 600);
129  CFinfo << "result[0] = " << result[0] << CFendl;
130  CFinfo << "result[1] = " << result[1] << CFendl << CFendl;
131 
132  test_coords = nodes.row(1);
133  ETYPE::compute_mapped_coordinate(test_coords, nodes, result);
134  vector_test(result,ETYPE::MappedCoordsT(1.0, -1.0),accumulator);
135  BOOST_CHECK_LT(boost::accumulators::max(accumulator.ulps), 600);
136  CFinfo << "result[0] = " << result[0] << CFendl;
137  CFinfo << "result[1] = " << result[1] << CFendl << CFendl;
138 
139  test_coords = nodes.row(2);
140  ETYPE::compute_mapped_coordinate(test_coords, nodes, result);
141  vector_test(result,ETYPE::MappedCoordsT(1.0, 1.0),accumulator);
142  BOOST_CHECK_LT(boost::accumulators::max(accumulator.ulps), 600);
143  CFinfo << "result[0] = " << result[0] << CFendl;
144  CFinfo << "result[1] = " << result[1] << CFendl << CFendl;
145 
146  test_coords = nodes.row(3);
147  ETYPE::compute_mapped_coordinate(test_coords, nodes, result);
148  vector_test(result,ETYPE::MappedCoordsT(-1.0, 1.0),accumulator);
149  BOOST_CHECK_LT(boost::accumulators::max(accumulator.ulps), 600);
150  CFinfo << "result[0] = " << result[0] << CFendl;
151  CFinfo << "result[1] = " << result[1] << CFendl << CFendl;
152 
153  // Try another element
154  ETYPE::NodesT nodes_2;
155  nodes_2 <<
156  1.0, 1.0,
157  2.0, 1.0,
158  2.0, 2.0,
159  1.0, 2.0;
160 
161  test_coords = nodes_2.row(0);
162  ETYPE::compute_mapped_coordinate(test_coords, nodes_2, result);
163  CFinfo << "result[0] = " << result[0] << CFendl;
164  CFinfo << "result[1] = " << result[1] << CFendl << CFendl;
165  vector_test(result,ETYPE::MappedCoordsT(-1.0, -1.0),accumulator);
166  BOOST_CHECK_LT(boost::accumulators::max(accumulator.ulps), 600);
167 
168  test_coords = nodes_2.row(1);
169  ETYPE::compute_mapped_coordinate(test_coords, nodes_2, result);
170  vector_test(result,ETYPE::MappedCoordsT(1.0, -1.0),accumulator);
171  BOOST_CHECK_LT(boost::accumulators::max(accumulator.ulps), 600);
172  CFinfo << "result[0] = " << result[0] << CFendl;
173  CFinfo << "result[1] = " << result[1] << CFendl << CFendl;
174 
175  test_coords = nodes_2.row(2);
176  ETYPE::compute_mapped_coordinate(test_coords, nodes_2, result);
177  vector_test(result,ETYPE::MappedCoordsT(1.0, 1.0),accumulator);
178  BOOST_CHECK_LT(boost::accumulators::max(accumulator.ulps), 600);
179  CFinfo << "result[0] = " << result[0] << CFendl;
180  CFinfo << "result[1] = " << result[1] << CFendl << CFendl;
181 
182  test_coords = nodes_2.row(3);
183  ETYPE::compute_mapped_coordinate(test_coords, nodes_2, result);
184  vector_test(result,ETYPE::MappedCoordsT(-1.0, 1.0),accumulator);
185  BOOST_CHECK_LT(boost::accumulators::max(accumulator.ulps), 600);
186  CFinfo << "result[0] = " << result[0] << CFendl;
187  CFinfo << "result[1] = " << result[1] << CFendl << CFendl;
188 
189 
190  // Try another element
191  const Real D=200./120.;
192  const Real eps = 0.0000000001;
193  ETYPE::NodesT nodes_3;
194  nodes_3 <<
195  -75. + eps, 0.,
196  -75.0+D, eps,
197  -75.0+D, D,
198  -75.0, D;
199 
200  test_coords << -75. + eps, 0.;
201  CFinfo << "test if " << test_coords.transpose() << " maps to (-1,-1) in element with nodes \n" << nodes_3 << CFendl;
202  ETYPE::compute_mapped_coordinate(test_coords, nodes_3, result);
203  CFinfo << "result[0] = " << result[0] << CFendl;
204  CFinfo << "result[1] = " << result[1] << CFendl << CFendl;
205  vector_test(result,ETYPE::MappedCoordsT(-1.0, -1.0),accumulator);
206  BOOST_CHECK_LT(boost::accumulators::max(accumulator.ulps), 300000);
207 
208 
209 }
210 
211 BOOST_AUTO_TEST_CASE( computeMappedGradient )
212 {
213  ETYPE::SF::GradientT expected;
214  const cf3::Real ksi = mapped_coords[0];
215  const cf3::Real eta = mapped_coords[1];
216  expected(0,0) = 0.25 * (-1 + eta);
217  expected(1,0) = 0.25 * (-1 + ksi);
218  expected(0,1) = 0.25 * ( 1 - eta);
219  expected(1,1) = 0.25 * (-1 - ksi);
220  expected(0,2) = 0.25 * ( 1 + eta);
221  expected(1,2) = 0.25 * ( 1 + ksi);
222  expected(0,3) = 0.25 * (-1 - eta);
223  expected(1,3) = 0.25 * ( 1 - ksi);
224  ETYPE::SF::GradientT result;
225  ETYPE::SF::compute_gradient(mapped_coords, result);
226  Accumulator accumulator;
227  vector_test(result, expected, accumulator);
228  BOOST_CHECK_LT(boost::accumulators::max(accumulator.ulps), 2);
229 }
230 
231 BOOST_AUTO_TEST_CASE( computeJacobianDeterminant )
232 {
233  // Shapefunction determinant at center should be a quarter of the cell volume
234  const ETYPE::CoordsT center_coords(0., 0.);
235  BOOST_CHECK_LT(boost::accumulators::max(test(4.*ETYPE::jacobian_determinant(center_coords, nodes), ETYPE::volume(nodes)).ulps), 1);
236 }
237 
238 BOOST_AUTO_TEST_CASE( computeJacobian )
239 {
240  ETYPE::JacobianT expected;
241  expected(0,0) = 0.2775;
242  expected(0,1) = -0.045;
243  expected(1,0) = 0.13625;
244  expected(1,1) = 0.5975;
245  ETYPE::JacobianT result;
246  ETYPE::compute_jacobian(mapped_coords, nodes, result);
247  Accumulator accumulator;
248  vector_test(result, expected, accumulator);
249  BOOST_CHECK_LT(boost::accumulators::max(accumulator.ulps), 15);
250 
251 
252  ETYPE::NodesT nodes;
253  nodes << 0. , 0. , 1. , 0., 1. , 1. , 0. , 1.;
254  ETYPE::MappedCoordsT d_mapped; d_mapped << 2. , 2.;
255  ETYPE::MappedCoordsT mapped; mapped << 0.5, 0.5;
256  ETYPE::JacobianT jacobian;
257  ETYPE::compute_jacobian(mapped, nodes, jacobian);
258 
259  RealVector2 d_phys = jacobian * d_mapped;
260  BOOST_CHECK_EQUAL( d_phys[0], 1. );
261  BOOST_CHECK_EQUAL( d_phys[1], 1. );
262 }
263 
264 BOOST_AUTO_TEST_CASE( computeJacobianAdjoint )
265 {
266  ETYPE::JacobianT expected;
267  expected(0,0) = 0.5975;
268  expected(0,1) = 0.045;
269  expected(1,0) = -0.13625;
270  expected(1,1) = 0.2775;
271  ETYPE::JacobianT result(2, 2);
272  ETYPE::compute_jacobian_adjoint(mapped_coords, nodes, result);
273  Accumulator accumulator;
274  vector_test(result, expected, accumulator);
275  BOOST_CHECK_LT(boost::accumulators::max(accumulator.ulps), 15);
276 }
277 
278 BOOST_AUTO_TEST_CASE( integrateConst )
279 {
280  // Shapefunction determinant should be double the volume for triangles
281  ConstFunctor ftor(nodes);
282  const Real vol = ETYPE::volume(nodes);
283 
284  cf3::Real result1 = 0.0;
285  cf3::Real result2 = 0.0;
286  cf3::Real result4 = 0.0;
287  cf3::Real result8 = 0.0;
288  cf3::Real result16 = 0.0;
289  cf3::Real result32 = 0.0;
290 
291  gauss_integrate<1, GeoShape::QUAD>(ftor, ftor.mapped_coords, result1);
292  gauss_integrate<2, GeoShape::QUAD>(ftor, ftor.mapped_coords, result2);
293  gauss_integrate<4, GeoShape::QUAD>(ftor, ftor.mapped_coords, result4);
294  gauss_integrate<8, GeoShape::QUAD>(ftor, ftor.mapped_coords, result8);
295  gauss_integrate<16, GeoShape::QUAD>(ftor, ftor.mapped_coords, result16);
296  gauss_integrate<32, GeoShape::QUAD>(ftor, ftor.mapped_coords, result32);
297 
298  BOOST_CHECK_LT(boost::accumulators::max(test(result1, vol).ulps), 1);
299  BOOST_CHECK_LT(boost::accumulators::max(test(result2, vol).ulps), 15);
300  BOOST_CHECK_LT(boost::accumulators::max(test(result4, vol).ulps), 15);
301  BOOST_CHECK_LT(boost::accumulators::max(test(result8, vol).ulps), 15);
302  BOOST_CHECK_LT(boost::accumulators::max(test(result16, vol).ulps), 15);
303  BOOST_CHECK_LT(boost::accumulators::max(test(result32, vol).ulps), 15);
304 }
305 
307 
308 BOOST_AUTO_TEST_SUITE_END()
309 
310 
cf3::solver::actions::Proto::jacobian
boost::proto::terminal< SFOp< JacobianOp > >::type const jacobian
Definition: ElementOperations.hpp:810
CFinfo
#define CFinfo
these are always defined
Definition: Log.hpp:104
cf3::mesh::LagrangeP1::Hexa3D::compute_jacobian_adjoint
static void compute_jacobian_adjoint(const MappedCoordsT &mapped_coord, const NodesT &nodes, JacobianT &result)
Definition: Hexa3D.cpp:308
ContinuousDictionary.hpp
rawmat::result2
cf3::Real result2[213]
Definition: utest-lss-test-matrix.hpp:656
BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(Volume)
Definition: utest-mesh-lagrangep1-quad2d.cpp:75
cf3::mesh::LagrangeP1::Hexa3D::compute_jacobian
static void compute_jacobian(const MappedCoordsT &mapped_coord, const NodesT &nodes, MatrixType &jacobian)
cf3::mesh::LagrangeP1::Hexa3D::jacobian_determinant
static Real jacobian_determinant(const MappedCoordsT &mapped_coord, const NodesT &nodes)
Definition: Hexa3D.cpp:169
Quad2D.hpp
Table.hpp
rawmat::result1
cf3::Real result1[201]
Definition: utest-lss-test-matrix.hpp:650
CFendl
#define CFendl
Definition: Log.hpp:109
cf3::solver::actions::Proto::max
Real max(const Real a, const Real b)
Maximum between two scalars.
Definition: Terminals.hpp:228
cf3::physics::Consts::e
Real e()
Definition of the Unit charge [C].
Definition: Consts.hpp:30
cf3::solver::actions::Proto::nodes
boost::proto::terminal< SFOp< NodesOp > >::type const nodes
Definition: ElementOperations.hpp:807
cf3::mesh::ElementTypeBase< Hexa3D, Hexa3D_traits >::NodesT
Eigen::Matrix< Real, nb_nodes, Hexa3D_traits::dimension > NodesT
Definition: ElementTypeBase.hpp:69
cf3::Tools::Testing::Accumulator::ulps
boost::accumulators::accumulator_set< Real, boost::accumulators::stats< boost::accumulators::tag::min, boost::accumulators::tag::mean, boost::accumulators::tag::max, boost::accumulators::tag::median, boost::accumulators::tag::lazy_variance > > ulps
Stores statistics for comparisons of inexact (floating-point) types using Units in the Last Place (UL...
Definition: Difference.hpp:55
Elements.hpp
cf3::mesh
Basic Classes for Mesh applications used by COOLFluiD.
Definition: AddPointRegion.cpp:33
cf3::mesh::LagrangeP1
namespace holding LagrangeP1 shape functions and elements
Definition: ElementTypes.hpp:28
cf3
Top-level namespace for coolfluid.
Definition: Action.cpp:18
Difference.hpp
Gauss.hpp
cf3::mesh::LagrangeP1::Hexa3D::compute_mapped_coordinate
static void compute_mapped_coordinate(const CoordsT &coord, const NodesT &nodes, MappedCoordsT &mapped_coord)
Definition: Hexa3D.cpp:85
cf3::mesh::ElementTypeBase< Hexa3D, Hexa3D_traits >::CoordsT
Eigen::Matrix< Real, Hexa3D_traits::dimension, 1 > CoordsT
Definition: ElementTypeBase.hpp:68
cf3::mesh::Integrators
Functions to provide integration over elements.
Definition: Gauss.hpp:26
LagrangeP1Quad2DFixture::mapped_coords
const ETYPE::MappedCoordsT mapped_coords
common values accessed by all tests goes here
Definition: utest-mesh-lagrangep1-quad2d.cpp:52
LagrangeP1Quad2DFixture::~LagrangeP1Quad2DFixture
~LagrangeP1Quad2DFixture()
common tear-down for each test case
Definition: utest-mesh-lagrangep1-quad2d.cpp:47
cf3::Uint
unsigned int Uint
typedef for unsigned int
Definition: CF.hpp:90
Log.hpp
cf3::mesh::LagrangeP1::Hexa3D::volume
static Real volume(const NodesT &nodes)
Definition: Hexa3D.cpp:324
cf3::RealVector2
Eigen::Matrix< Real, 2, 1 > RealVector2
Fixed size 2x1 column vector.
Definition: MatrixTypes.hpp:40
cf3::Tools::Testing::vector_test
void vector_test(const VectorT &A, const VectorT &B, Accumulator &Result)
Compares vector-like sequences.
Definition: Difference.hpp:139
LagrangeP1Quad2DFixture::LagrangeP1Quad2DFixture
LagrangeP1Quad2DFixture()
common setup for each test case
Definition: utest-mesh-lagrangep1-quad2d.cpp:39
cf3::mesh::ElementTypeBase< Hexa3D, Hexa3D_traits >::JacobianT
Eigen::Matrix< Real, Hexa3D_traits::SF::dimensionality, Hexa3D_traits::dimension > JacobianT
Definition: ElementTypeBase.hpp:70
cf3::common
Most basic kernel library.
Definition: Action.cpp:19
cf3::mesh::LagrangeP1::Quad2D
2D Lagrange P1 Quadrilateral Element type This class provides the lagrangian shape function describin...
Definition: Quad2D.hpp:37
cf3::Tools::Testing::Accumulator
Stores the results of the difference::test() function.
Definition: Difference.hpp:49
Generated on Sun Jun 14 2015 21:20:17 for COOLFluiD by doxygen 1.8.9.1
Send comments to:
COOLFluiD Web Admin