Reference documentation for deal.II version Git e3a3ec7800 2020-08-07 14:08:19 +0200
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
grid_tools_nontemplates.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/point.h>
17 
19 
20 #include <vector>
21 
22 // GridTools functions that are template specializations (i.e., only compiled
23 // once without expand_instantiations)
24 
26 
27 
28 namespace GridTools
29 {
30  template <>
31  double
32  cell_measure<1>(const std::vector<Point<1>> & all_vertices,
34  {
36 
37  return all_vertices[vertex_indices[1]][0] -
38  all_vertices[vertex_indices[0]][0];
39  }
40 
41 
42 
43  template <>
44  double
45  cell_measure<2>(const std::vector<Point<2>> & all_vertices,
47  {
48  if (vertex_indices.size() == 3) // triangle
49  {
50  const double x[3] = {all_vertices[vertex_indices[0]](0),
51  all_vertices[vertex_indices[1]](0),
52  all_vertices[vertex_indices[2]](0)};
53 
54  const double y[3] = {all_vertices[vertex_indices[0]](1),
55  all_vertices[vertex_indices[1]](1),
56  all_vertices[vertex_indices[2]](1)};
57 
58  return 0.5 * std::abs((x[0] - x[2]) * (y[1] - y[0]) -
59  (x[0] - x[1]) * (y[2] - y[0]));
60  }
61 
63 
64  /*
65  Get the computation of the measure by this little Maple script. We
66  use the blinear mapping of the unit quad to the real quad. However,
67  every transformation mapping the unit faces to straight lines should
68  do.
69 
70  Remember that the area of the quad is given by
71  \int_K 1 dx dy = \int_{\hat K} |det J| d(xi) d(eta)
72 
73  # x and y are arrays holding the x- and y-values of the four vertices
74  # of this cell in real space.
75  x := array(0..3);
76  y := array(0..3);
77  z := array(0..3);
78  tphi[0] := (1-xi)*(1-eta):
79  tphi[1] := xi*(1-eta):
80  tphi[2] := (1-xi)*eta:
81  tphi[3] := xi*eta:
82  x_real := sum(x[s]*tphi[s], s=0..3):
83  y_real := sum(y[s]*tphi[s], s=0..3):
84  z_real := sum(z[s]*tphi[s], s=0..3):
85 
86  Jxi := <diff(x_real,xi) | diff(y_real,xi) | diff(z_real,xi)>;
87  Jeta := <diff(x_real,eta)| diff(y_real,eta)| diff(z_real,eta)>;
88  with(VectorCalculus):
89  J := CrossProduct(Jxi, Jeta);
90  detJ := sqrt(J[1]^2 + J[2]^2 +J[3]^2);
91 
92  # measure := evalf (Int (Int (detJ, xi=0..1, method = _NCrule ) ,
93  eta=0..1, method = _NCrule ) ): # readlib(C):
94 
95  # C(measure, optimized);
96 
97  additional optimizaton: divide by 2 only one time
98  */
99 
100  const double x[4] = {all_vertices[vertex_indices[0]](0),
101  all_vertices[vertex_indices[1]](0),
102  all_vertices[vertex_indices[2]](0),
103  all_vertices[vertex_indices[3]](0)};
104 
105  const double y[4] = {all_vertices[vertex_indices[0]](1),
106  all_vertices[vertex_indices[1]](1),
107  all_vertices[vertex_indices[2]](1),
108  all_vertices[vertex_indices[3]](1)};
109 
110  return (-x[1] * y[0] + x[1] * y[3] + y[0] * x[2] + x[0] * y[1] -
111  x[0] * y[2] - y[1] * x[3] - x[2] * y[3] + x[3] * y[2]) /
112  2;
113  }
114 
115 
116 
117  template <>
118  double
119  cell_measure<3>(const std::vector<Point<3>> & all_vertices,
121  {
122  if (vertex_indices.size() == 4) // tetrahedron
123  {
124  const auto &a = all_vertices[vertex_indices[0]];
125  const auto &b = all_vertices[vertex_indices[1]];
126  const auto &c = all_vertices[vertex_indices[2]];
127  const auto &d = all_vertices[vertex_indices[3]];
128 
129  return (1.0 / 6.0) * std::abs((a - d) * cross_product_3d(b - d, c - d));
130  }
131 
133  // note that this is the
134  // cell_measure based on the new
135  // deal.II numbering. When called
136  // from inside GridReordering make
137  // sure that you reorder the
138  // vertex_indices before
139  const double x[8] = {all_vertices[vertex_indices[0]](0),
140  all_vertices[vertex_indices[1]](0),
141  all_vertices[vertex_indices[2]](0),
142  all_vertices[vertex_indices[3]](0),
143  all_vertices[vertex_indices[4]](0),
144  all_vertices[vertex_indices[5]](0),
145  all_vertices[vertex_indices[6]](0),
146  all_vertices[vertex_indices[7]](0)};
147  const double y[8] = {all_vertices[vertex_indices[0]](1),
148  all_vertices[vertex_indices[1]](1),
149  all_vertices[vertex_indices[2]](1),
150  all_vertices[vertex_indices[3]](1),
151  all_vertices[vertex_indices[4]](1),
152  all_vertices[vertex_indices[5]](1),
153  all_vertices[vertex_indices[6]](1),
154  all_vertices[vertex_indices[7]](1)};
155  const double z[8] = {all_vertices[vertex_indices[0]](2),
156  all_vertices[vertex_indices[1]](2),
157  all_vertices[vertex_indices[2]](2),
158  all_vertices[vertex_indices[3]](2),
159  all_vertices[vertex_indices[4]](2),
160  all_vertices[vertex_indices[5]](2),
161  all_vertices[vertex_indices[6]](2),
162  all_vertices[vertex_indices[7]](2)};
163 
164  /*
165  This is the same Maple script as in the barycenter method above
166  except of that here the shape functions tphi[0]-tphi[7] are ordered
167  according to the lexicographic numbering.
168 
169  x := array(0..7):
170  y := array(0..7):
171  z := array(0..7):
172  tphi[0] := (1-xi)*(1-eta)*(1-zeta):
173  tphi[1] := xi*(1-eta)*(1-zeta):
174  tphi[2] := (1-xi)* eta*(1-zeta):
175  tphi[3] := xi* eta*(1-zeta):
176  tphi[4] := (1-xi)*(1-eta)*zeta:
177  tphi[5] := xi*(1-eta)*zeta:
178  tphi[6] := (1-xi)* eta*zeta:
179  tphi[7] := xi* eta*zeta:
180  x_real := sum(x[s]*tphi[s], s=0..7):
181  y_real := sum(y[s]*tphi[s], s=0..7):
182  z_real := sum(z[s]*tphi[s], s=0..7):
183  with (linalg):
184  J := matrix(3,3, [[diff(x_real, xi), diff(x_real, eta), diff(x_real,
185  zeta)], [diff(y_real, xi), diff(y_real, eta), diff(y_real, zeta)],
186  [diff(z_real, xi), diff(z_real, eta), diff(z_real, zeta)]]):
187  detJ := det (J):
188 
189  measure := simplify ( int ( int ( int (detJ, xi=0..1), eta=0..1),
190  zeta=0..1)):
191 
192  readlib(C):
193 
194  C(measure, optimized);
195 
196  The C code produced by this maple script is further optimized by
197  hand. In particular, division by 12 is performed only once, not
198  hundred of times.
199  */
200 
201  const double t3 = y[3] * x[2];
202  const double t5 = z[1] * x[5];
203  const double t9 = z[3] * x[2];
204  const double t11 = x[1] * y[0];
205  const double t14 = x[4] * y[0];
206  const double t18 = x[5] * y[7];
207  const double t20 = y[1] * x[3];
208  const double t22 = y[5] * x[4];
209  const double t26 = z[7] * x[6];
210  const double t28 = x[0] * y[4];
211  const double t34 =
212  z[3] * x[1] * y[2] + t3 * z[1] - t5 * y[7] + y[7] * x[4] * z[6] +
213  t9 * y[6] - t11 * z[4] - t5 * y[3] - t14 * z[2] + z[1] * x[4] * y[0] -
214  t18 * z[3] + t20 * z[0] - t22 * z[0] - y[0] * x[5] * z[4] - t26 * y[3] +
215  t28 * z[2] - t9 * y[1] - y[1] * x[4] * z[0] - t11 * z[5];
216  const double t37 = y[1] * x[0];
217  const double t44 = x[1] * y[5];
218  const double t46 = z[1] * x[0];
219  const double t49 = x[0] * y[2];
220  const double t52 = y[5] * x[7];
221  const double t54 = x[3] * y[7];
222  const double t56 = x[2] * z[0];
223  const double t58 = x[3] * y[2];
224  const double t64 = -x[6] * y[4] * z[2] - t37 * z[2] + t18 * z[6] -
225  x[3] * y[6] * z[2] + t11 * z[2] + t5 * y[0] +
226  t44 * z[4] - t46 * y[4] - t20 * z[7] - t49 * z[6] -
227  t22 * z[1] + t52 * z[3] - t54 * z[2] - t56 * y[4] -
228  t58 * z[0] + y[1] * x[2] * z[0] + t9 * y[7] + t37 * z[4];
229  const double t66 = x[1] * y[7];
230  const double t68 = y[0] * x[6];
231  const double t70 = x[7] * y[6];
232  const double t73 = z[5] * x[4];
233  const double t76 = x[6] * y[7];
234  const double t90 = x[4] * z[0];
235  const double t92 = x[1] * y[3];
236  const double t95 = -t66 * z[3] - t68 * z[2] - t70 * z[2] + t26 * y[5] -
237  t73 * y[6] - t14 * z[6] + t76 * z[2] - t3 * z[6] +
238  x[6] * y[2] * z[4] - z[3] * x[6] * y[2] + t26 * y[4] -
239  t44 * z[3] - x[1] * y[2] * z[0] + x[5] * y[6] * z[4] +
240  t54 * z[5] + t90 * y[2] - t92 * z[2] + t46 * y[2];
241  const double t102 = x[2] * y[0];
242  const double t107 = y[3] * x[7];
243  const double t114 = x[0] * y[6];
244  const double t125 =
245  y[0] * x[3] * z[2] - z[7] * x[5] * y[6] - x[2] * y[6] * z[4] +
246  t102 * z[6] - t52 * z[6] + x[2] * y[4] * z[6] - t107 * z[5] - t54 * z[6] +
247  t58 * z[6] - x[7] * y[4] * z[6] + t37 * z[5] - t114 * z[4] + t102 * z[4] -
248  z[1] * x[2] * y[0] + t28 * z[6] - y[5] * x[6] * z[4] -
249  z[5] * x[1] * y[4] - t73 * y[7];
250  const double t129 = z[0] * x[6];
251  const double t133 = y[1] * x[7];
252  const double t145 = y[1] * x[5];
253  const double t156 = t90 * y[6] - t129 * y[4] + z[7] * x[2] * y[6] -
254  t133 * z[5] + x[5] * y[3] * z[7] - t26 * y[2] -
255  t70 * z[3] + t46 * y[3] + z[5] * x[7] * y[4] +
256  z[7] * x[3] * y[6] - t49 * z[4] + t145 * z[7] -
257  x[2] * y[7] * z[6] + t70 * z[5] + t66 * z[5] -
258  z[7] * x[4] * y[6] + t18 * z[4] + x[1] * y[4] * z[0];
259  const double t160 = x[5] * y[4];
260  const double t165 = z[1] * x[7];
261  const double t178 = z[1] * x[3];
262  const double t181 =
263  t107 * z[6] + t22 * z[7] + t76 * z[3] + t160 * z[1] - x[4] * y[2] * z[6] +
264  t70 * z[4] + t165 * y[5] + x[7] * y[2] * z[6] - t76 * z[5] - t76 * z[4] +
265  t133 * z[3] - t58 * z[1] + y[5] * x[0] * z[4] + t114 * z[2] - t3 * z[7] +
266  t20 * z[2] + t178 * y[7] + t129 * y[2];
267  const double t207 = t92 * z[7] + t22 * z[6] + z[3] * x[0] * y[2] -
268  x[0] * y[3] * z[2] - z[3] * x[7] * y[2] - t165 * y[3] -
269  t9 * y[0] + t58 * z[7] + y[3] * x[6] * z[2] +
270  t107 * z[2] + t73 * y[0] - x[3] * y[5] * z[7] +
271  t3 * z[0] - t56 * y[6] - z[5] * x[0] * y[4] +
272  t73 * y[1] - t160 * z[6] + t160 * z[0];
273  const double t228 = -t44 * z[7] + z[5] * x[6] * y[4] - t52 * z[4] -
274  t145 * z[4] + t68 * z[4] + t92 * z[5] - t92 * z[0] +
275  t11 * z[3] + t44 * z[0] + t178 * y[5] - t46 * y[5] -
276  t178 * y[0] - t145 * z[0] - t20 * z[5] - t37 * z[3] -
277  t160 * z[7] + t145 * z[3] + x[4] * y[6] * z[2];
278 
279  return (t34 + t64 + t95 + t125 + t156 + t181 + t207 + t228) / 12.;
280  }
281 
282 
283 
284  namespace
285  {
286  // the following class is only
287  // needed in 2d, so avoid trouble
288  // with compilers warning otherwise
289  class Rotate2d
290  {
291  public:
292  explicit Rotate2d(const double angle)
293  : angle(angle)
294  {}
295  Point<2>
296  operator()(const Point<2> &p) const
297  {
298  return {std::cos(angle) * p(0) - std::sin(angle) * p(1),
299  std::sin(angle) * p(0) + std::cos(angle) * p(1)};
300  }
301 
302  private:
303  const double angle;
304  };
305  } // namespace
306 
307 
308 
309  template <>
310  void
312  {
313  transform(Rotate2d(angle), triangulation);
314  }
315 
316 
317 
318  template <>
319  void
321  {
322  (void)angle;
323  (void)triangulation;
324 
325  AssertThrow(
326  false, ExcMessage("GridTools::rotate() is not available for dim = 3."));
327  }
328 } /* namespace GridTools */
329 
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1568
std::vector< unsigned int > vertex_indices
Definition: tria.cc:2244
#define AssertThrow(cond, exc)
Definition: exceptions.h:1521
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2430
const double angle
static ::ExceptionBase & ExcMessage(std::string arg1)
double cell_measure< 3 >(const std::vector< Point< 3 >> &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void rotate(const double angle, Triangulation< dim > &triangulation)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
double cell_measure< 1 >(const std::vector< Point< 1 >> &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
double cell_measure< 2 >(const std::vector< Point< 2 >> &all_vertices, const ArrayView< const unsigned int > &vertex_indices)