Flow123d  master-f44eb46
element_values.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2019 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file fe_values.cc
15  * @brief Class ElementValues calculates data related to transformation
16  * of reference cell to actual cell (Jacobian, inverse Jacobian,
17  * determinant, point coordinates etc.).
18  * @author Jan Stebel
19  */
20 
21 #include "fem/mapping_p1.hh"
22 #include "quadrature/quadrature.hh"
23 #include "fem/element_values.hh"
24 
25 
26 
27 using namespace arma;
28 using namespace std;
29 
30 
31 
32 
33 
34 
35 
36 
38  : n_points(np)
39 {
40  bar_coords.resize(np);
41  weights.resize(np);
42 }
43 
44 
45 template<unsigned int spacedim>
46 ElementData<spacedim>::ElementData(unsigned int size,
47  UpdateFlags flags,
48  unsigned int dim)
49 : dim_(dim),
50  JxW_values(flags & update_JxW_values ? size : 0),
51  side_JxW_values(flags & update_side_JxW_values ? size : 0),
52  jacobians(spacedim, dim, flags & update_jacobians ? size : 0),
53  determinants(flags & update_volume_elements ? size : 0),
54  inverse_jacobians(dim, spacedim, flags & update_inverse_jacobians ? size : 0),
55  points(spacedim, 1, flags & update_quadrature_points ? size : 0),
56  normal_vectors(spacedim, 1, flags & update_normal_vectors ? size : 0),
57  update_flags(flags)
58 {}
59 
60 
61 template<unsigned int spacedim>
63 {
64  if (cell.is_valid() || side.is_valid())
65  {
66  if (cell.is_valid())
67  printf("cell %d dim %d ", cell.elm_idx(), cell.dim());
68  else if (side.is_valid())
69  printf("cell %d dim %d side %d ", side.elem_idx(), side.dim(), side.side_idx());
70 
71  printf(" det[");
72  for (auto d : determinants) printf("%g ", d); printf("]");
73 
74  printf(" JxW[");
75  for (auto j : JxW_values) printf("%g ", j); printf("]");
76 
77  printf(" side_JxW[");
78  for (auto j : side_JxW_values) printf("%g ", j); printf("]");
79 
80  printf(" nv[");
81  for (unsigned int i=0; i<normal_vectors.size(); i++)
82  {
83  auto n = normal_vectors.vec<spacedim>(i);
84  printf(" [");
85  for (unsigned int c=0; c<spacedim; c++) printf("%g ", n[c]);
86  printf("]");
87  }
88  printf("]\n");
89  }
90 }
91 
92 
93 
94 
95 
96 template<unsigned int spacedim>
98 {
99  ASSERT( q.dim() == dim_ );
100  ASSERT( q.size() == n_points_ );
101  RefElementData *ref_data = new RefElementData(q.size());
102 
103  for (unsigned int i=0; i<q.size(); i++)
104  {
105  switch (q.dim())
106  {
107  case 0:
108  ref_data->bar_coords[i] = arma::vec( "0" );
109  //ref_data->bar_coords[i] = RefElement<0>::local_to_bary(q.point<0>(i));
110  break;
111  case 1:
112  ref_data->bar_coords[i] = RefElement<1>::local_to_bary(q.point<1>(i));
113  break;
114  case 2:
115  ref_data->bar_coords[i] = RefElement<2>::local_to_bary(q.point<2>(i));
116  break;
117  case 3:
118  ref_data->bar_coords[i] = RefElement<3>::local_to_bary(q.point<3>(i));
119  break;
120  default:
121  ASSERT_PERMANENT(false)(q.dim()).error("Unsupported dimension.\n");
122  break;
123  }
124  ref_data->weights[i] = q.weight(i);
125  }
126 
127  return ref_data;
128 }
129 
130 
131 
132 template<unsigned int spacedim>
134 {
135  switch (dim_)
136  {
137  case 0:
138  flags = MappingP1<0,spacedim>::update_each(flags);
139  break;
140  case 1:
141  flags = MappingP1<1,spacedim>::update_each(flags);
142  break;
143  case 2:
144  flags = MappingP1<2,spacedim>::update_each(flags);
145  break;
146  case 3:
147  flags = MappingP1<3,spacedim>::update_each(flags);
148  break;
149  default:
150  ASSERT_PERMANENT(false)(dim_).error("Unsupported dimension.\n");
151  break;
152  }
153  return flags;
154 }
155 
156 
157 
158 template<unsigned int spacedim>
160  Quadrature &_quadrature,
161  UpdateFlags _flags,
162  unsigned int dim)
163 : dim_(dim),
164  n_points_(_quadrature.size()),
165  n_sides_(_quadrature.dim() == dim ? 0 : dim+1),
166  ref_data(nullptr),
167  side_ref_data(n_sides_),
168  data(n_points_, update_each(_flags), dim)
169 {
170  //if (dim == 0) return; // avoid unnecessary allocation of dummy 0 dimensional objects
171  if ( _quadrature.dim() == dim )
172  {
173  // precompute element data
174  ref_data = init_ref_data(_quadrature);
175  }
176  else if ( _quadrature.dim() + 1 == dim )
177  {
178  // precompute side data
179  for (unsigned int sid = 0; sid < n_sides_; sid++)
180  {
181  Quadrature side_quad(dim);
182  // transform the side quadrature points to the cell quadrature points
183  switch (dim)
184  {
185  case 1:
186  side_quad = _quadrature.make_from_side<1>(sid);
187  break;
188  case 2:
189  side_quad = _quadrature.make_from_side<2>(sid);
190  break;
191  case 3:
192  side_quad = _quadrature.make_from_side<3>(sid);
193  break;
194  default:
195  ASSERT_PERMANENT(false)(dim).error("Unsupported dimension.\n");
196  break;
197  }
198  side_ref_data[sid] = init_ref_data(side_quad);
199  }
200  }
201 }
202 
203 
204 template<unsigned int spacedim>
206 {
207  if (ref_data) delete ref_data;
208 
209  for (unsigned int sid=0; sid<n_sides_; sid++)
210  delete side_ref_data[sid];
211 }
212 
213 
214 template<unsigned int spacedim>
216 {
217  ASSERT_EQ( dim_, cell.dim() );
218  data.cell = cell;
219 
220  // calculate Jacobian of mapping, JxW, inverse Jacobian
221  switch (dim_)
222  {
223  case 0:
224  if (cell.is_valid() && data.update_flags & update_quadrature_points)
225  data.points.set(0) = Armor::vec<spacedim>( MappingP1<0, spacedim>::element_map(cell) );
226  break;
227  case 1:
228  fill_data<1>();
229  break;
230  case 2:
231  fill_data<2>();
232  break;
233  case 3:
234  fill_data<3>();
235  break;
236  default:
237  ASSERT_PERMANENT(false)(dim_).error("Unsupported dimension.\n");
238  break;
239  }
240 }
241 
242 
243 template<unsigned int spacedim>
244 void ElementValues<spacedim>::reinit(const Side & cell_side)
245 {
246  ASSERT_EQ( dim_, cell_side.dim()+1 );
247  data.side = cell_side;
248 
249  // calculate Jacobian of mapping, JxW, inverse Jacobian, normal vector(s)
250  switch (dim_)
251  {
252  case 1:
253  fill_data<1>();
254  fill_side_data<1>();
255  break;
256  case 2:
257  fill_data<2>();
258  fill_side_data<2>();
259  break;
260  case 3:
261  fill_data<3>();
262  fill_side_data<3>();
263  break;
264  default:
265  ASSERT_PERMANENT(false)(dim_).error("Unsupported dimension.\n");
266  break;
267  }
268 }
269 
270 
271 template<unsigned int spacedim>
272 template<unsigned int dim>
274 {
275  typename MappingP1<dim,spacedim>::ElementMap coords;
276 
277  if ((data.update_flags & update_jacobians) |
278  (data.update_flags & update_volume_elements) |
279  (data.update_flags & update_JxW_values) |
280  (data.update_flags & update_inverse_jacobians) |
281  (data.update_flags & update_normal_vectors) |
282  (data.update_flags & update_quadrature_points))
283  {
284  if (cell().is_valid())
285  coords = MappingP1<dim,spacedim>::element_map(cell());
286  else
287  coords = MappingP1<dim,spacedim>::element_map(side().element());
288  }
289 
290  // calculation of Jacobian dependent data
291  if ((data.update_flags & update_jacobians) |
292  (data.update_flags & update_volume_elements) |
293  (data.update_flags & update_JxW_values) |
294  (data.update_flags & update_inverse_jacobians) |
295  (data.update_flags & update_normal_vectors))
296  {
297  arma::mat::fixed<spacedim,dim> jac = MappingP1<dim,spacedim>::jacobian(coords);
298 
299  // update Jacobians
300  if (data.update_flags & update_jacobians)
301  for (unsigned int i=0; i<n_points_; i++)
302  data.jacobians.set(i) = Armor::mat<spacedim,dim>( jac );
303 
304  // calculation of determinant dependent data
305  if ((data.update_flags & update_volume_elements) |
306  (data.update_flags & update_JxW_values))
307  {
308  double det = fabs(::determinant(jac));
309 
310  // update determinants
311  if (data.update_flags & update_volume_elements)
312  for (unsigned int i=0; i<n_points_; i++)
313  data.determinants[i] = det;
314 
315  // update JxW values
316  if (data.update_flags & update_JxW_values)
317  for (unsigned int i=0; i<n_points_; i++)
318  data.JxW_values[i] = det*ref_data->weights[i];
319  }
320 
321  // update inverse Jacobians
322  if (data.update_flags & update_inverse_jacobians)
323  {
324  arma::mat::fixed<dim,spacedim> ijac;
325  if (dim==spacedim)
326  {
327  ijac = inv(jac);
328  }
329  else
330  {
331  ijac = pinv(jac);
332  }
333  for (unsigned int i=0; i<n_points_; i++)
334  data.inverse_jacobians.set(i) = Armor::mat<dim,spacedim>( ijac );
335  }
336  }
337 
338  // quadrature points in the actual cell coordinate system
339  if (cell().is_valid() && data.update_flags & update_quadrature_points)
340  {
341  for (unsigned int i=0; i<n_points_; i++)
342  data.points.set(i) = Armor::vec<spacedim>( coords*ref_data->bar_coords[i] );
343  }
344 }
345 
346 
347 template<unsigned int spacedim>
348 template<unsigned int dim>
350 {
351  const unsigned int side_idx = side().side_idx();
352 
353  // calculation of normal vectors to the side
354  if (data.update_flags & update_normal_vectors)
355  {
356  arma::vec::fixed<spacedim> n_cell;
357  n_cell = trans(data.inverse_jacobians.template mat<dim,spacedim>(0))*RefElement<dim>::normal_vector(side_idx);
358  n_cell = n_cell/norm(n_cell,2);
359  for (unsigned int i=0; i<n_points_; i++)
360  data.normal_vectors.set(i) = Armor::vec<spacedim>( n_cell );
361  }
362 
363  // Quadrature points in the actual cell coordinate system.
364  // The points location can vary from side to side.
365  if (data.update_flags & update_quadrature_points)
366  {
368  for (unsigned int i=0; i<n_points_; i++)
369  data.points.set(i) = Armor::vec<spacedim>( coords*side_ref_data[side_idx]->bar_coords[i] );
370  }
371 
372  if (data.update_flags & update_side_JxW_values)
373  {
374  double side_det;
375  if (dim <= 1)
376  {
377  side_det = 1;
378  }
379  else
380  {
381  arma::mat::fixed<spacedim,dim> side_coords;
382  arma::mat::fixed<spacedim, MatrixSizes<dim>::dim_minus_one > side_jac; // some compilers complain for case dim==0
383 
384  // calculation of side Jacobian
385  for (unsigned int n=0; n<dim; n++)
386  for (unsigned int c=0; c<spacedim; c++)
387  side_coords(c,n) = (*data.side.node(n))[c];
388  side_jac = MappingP1<MatrixSizes<dim>::dim_minus_one,spacedim>::jacobian(side_coords);
389 
390  // calculation of JxW
391  side_det = fabs(::determinant(side_jac));
392  }
393  for (unsigned int i=0; i<n_points_; i++)
394  data.side_JxW_values[i] = side_det*side_ref_data[side_idx]->weights[i];
395  }
396 }
397 
398 
399 
400 
401 
402 template class ElementValues<3>;
403 
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:348
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
ArmaVec< Type, nr > vec(uint mat_index) const
Definition: armor.hh:821
unsigned int size() const
Definition: armor.hh:728
bool is_valid() const
Definition: accessors.hh:184
unsigned int dim() const
Definition: accessors.hh:188
std::vector< double > determinants
Determinants of Jacobians at quadrature points.
ElementAccessor< spacedim > cell
Iterator to last updated cell.
ElementData(unsigned int size, UpdateFlags flags, unsigned int dim)
Resize the data arrays.
void print()
Print calculated data.
std::vector< double > JxW_values
Transformed quadrature weights.
Armor::array normal_vectors
Normal vectors (spacedim) to the element at the quadrature points lying on a side.
std::vector< double > side_JxW_values
JxW values for sides.
Side side
Iterator to last updated cell side.
Class for computation of data on cell and side.
const unsigned int n_sides_
Number of sides in reference cell.
RefElementData * ref_data
Data on reference element.
ElementValues(Quadrature &_quadrature, UpdateFlags _flags, unsigned int dim)
Constructor.
RefElementData * init_ref_data(const Quadrature &q)
Precompute data on reference element.
void fill_side_data()
Calculates the mapping data on a side of a cell.
~ElementValues()
Correct deallocation of objects created by 'initialize' methods.
std::vector< RefElementData * > side_ref_data
Data on reference element (for each side ).
void reinit(const ElementAccessor< spacedim > &cell)
Update cell-dependent data (gradients, Jacobians etc.)
UpdateFlags update_each(UpdateFlags flags)
Determine quantities to be recomputed on each cell.
void fill_data()
Compute data from reference cell and using MappingP1.
Affine mapping between reference and actual cell.
Definition: mapping_p1.hh:59
static arma::mat::fixed< spacedim, dim > jacobian(const ElementMap &coords)
Definition: mapping_p1.cc:58
arma::mat::fixed< spacedim, dim+1 > ElementMap
Definition: mapping_p1.hh:64
static ElementMap element_map(ElementAccessor< 3 > elm)
Definition: mapping_p1.cc:48
static UpdateFlags update_each(UpdateFlags flags)
Determines which additional quantities have to be computed.
Definition: mapping_p1.cc:30
Base class for quadrature rules on simplices in arbitrary dimensions.
Definition: quadrature.hh:48
Quadrature make_from_side(unsigned int sid) const
Definition: quadrature.cc:47
unsigned int size() const
Returns number of quadrature points.
Definition: quadrature.hh:86
unsigned int dim() const
Definition: quadrature.hh:72
Armor::ArmaVec< double, point_dim > point(unsigned int i) const
Returns the ith quadrature point.
Definition: quadrature.hh:151
double weight(unsigned int i) const
Returns the ith weight.
Definition: quadrature.hh:103
Structure for storing the precomputed element data.
RefElementData(unsigned int np)
Resize vectors to size np.
std::vector< arma::vec > bar_coords
Barycentric coordinates of quadrature points.
std::vector< double > weights
Quadrature weights.
static BaryPoint local_to_bary(const LocalPoint &lp)
Converts from local to barycentric coordinates.
Definition: ref_element.cc:187
static LocalPoint normal_vector(unsigned int sid)
unsigned int elem_idx() const
Returns index of element in Mesh::element_vec_.
Definition: accessors.hh:448
unsigned int side_idx() const
Returns local index of the side on the element.
Definition: accessors.hh:444
bool is_valid() const
Returns true if the side has assigned element.
Definition: accessors.hh:452
unsigned int dim() const
Returns dimension of the side, that is dimension of the element minus one.
Class ElementValues calculates data related to transformation of reference cell to actual cell (Jacob...
double determinant(const T &M)
Calculates determinant of a rectangular matrix.
Class MappingP1 implements the affine transformation of the unit cell onto the actual cell.
ArmaMat< double, N, M > mat
Definition: armor.hh:888
ArmaVec< double, N > vec
Definition: armor.hh:885
void printf(BasicWriter< Char > &w, BasicCStringRef< Char > format, ArgList args)
Definition: printf.h:444
Basic definitions of numerical quadrature rules.
UpdateFlags
Enum type UpdateFlags indicates which quantities are to be recomputed on each finite element cell.
Definition: update_flags.hh:68
@ update_volume_elements
Determinant of the Jacobian.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_jacobians
Volume element.
@ update_side_JxW_values
Transformed quadrature weight for cell sides.
@ update_inverse_jacobians
Volume element.
@ update_quadrature_points
Transformed quadrature points.