Reference documentation for deal.II version 9.2.0
\(\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
33  const std::vector<Point<1>> &all_vertices,
35  {
36  return all_vertices[vertex_indices[1]][0] -
37  all_vertices[vertex_indices[0]][0];
38  }
39 
40 
41 
42  template <>
43  double
45  const std::vector<Point<2>> &all_vertices,
47  {
48  /*
49  Get the computation of the measure by this little Maple script. We
50  use the blinear mapping of the unit quad to the real quad. However,
51  every transformation mapping the unit faces to straight lines should
52  do.
53 
54  Remember that the area of the quad is given by
55  \int_K 1 dx dy = \int_{\hat K} |det J| d(xi) d(eta)
56 
57  # x and y are arrays holding the x- and y-values of the four vertices
58  # of this cell in real space.
59  x := array(0..3);
60  y := array(0..3);
61  z := array(0..3);
62  tphi[0] := (1-xi)*(1-eta):
63  tphi[1] := xi*(1-eta):
64  tphi[2] := (1-xi)*eta:
65  tphi[3] := xi*eta:
66  x_real := sum(x[s]*tphi[s], s=0..3):
67  y_real := sum(y[s]*tphi[s], s=0..3):
68  z_real := sum(z[s]*tphi[s], s=0..3):
69 
70  Jxi := <diff(x_real,xi) | diff(y_real,xi) | diff(z_real,xi)>;
71  Jeta := <diff(x_real,eta)| diff(y_real,eta)| diff(z_real,eta)>;
72  with(VectorCalculus):
73  J := CrossProduct(Jxi, Jeta);
74  detJ := sqrt(J[1]^2 + J[2]^2 +J[3]^2);
75 
76  # measure := evalf (Int (Int (detJ, xi=0..1, method = _NCrule ) ,
77  eta=0..1, method = _NCrule ) ): # readlib(C):
78 
79  # C(measure, optimized);
80 
81  additional optimizaton: divide by 2 only one time
82  */
83 
84  const double x[4] = {all_vertices[vertex_indices[0]](0),
85  all_vertices[vertex_indices[1]](0),
86  all_vertices[vertex_indices[2]](0),
87  all_vertices[vertex_indices[3]](0)};
88 
89  const double y[4] = {all_vertices[vertex_indices[0]](1),
90  all_vertices[vertex_indices[1]](1),
91  all_vertices[vertex_indices[2]](1),
92  all_vertices[vertex_indices[3]](1)};
93 
94  return (-x[1] * y[0] + x[1] * y[3] + y[0] * x[2] + x[0] * y[1] -
95  x[0] * y[2] - y[1] * x[3] - x[2] * y[3] + x[3] * y[2]) /
96  2;
97  }
98 
99 
100 
101  template <>
102  double
104  const std::vector<Point<3>> &all_vertices,
106  {
107  // note that this is the
108  // cell_measure based on the new
109  // deal.II numbering. When called
110  // from inside GridReordering make
111  // sure that you reorder the
112  // vertex_indices before
113  const double x[8] = {all_vertices[vertex_indices[0]](0),
114  all_vertices[vertex_indices[1]](0),
115  all_vertices[vertex_indices[2]](0),
116  all_vertices[vertex_indices[3]](0),
117  all_vertices[vertex_indices[4]](0),
118  all_vertices[vertex_indices[5]](0),
119  all_vertices[vertex_indices[6]](0),
120  all_vertices[vertex_indices[7]](0)};
121  const double y[8] = {all_vertices[vertex_indices[0]](1),
122  all_vertices[vertex_indices[1]](1),
123  all_vertices[vertex_indices[2]](1),
124  all_vertices[vertex_indices[3]](1),
125  all_vertices[vertex_indices[4]](1),
126  all_vertices[vertex_indices[5]](1),
127  all_vertices[vertex_indices[6]](1),
128  all_vertices[vertex_indices[7]](1)};
129  const double z[8] = {all_vertices[vertex_indices[0]](2),
130  all_vertices[vertex_indices[1]](2),
131  all_vertices[vertex_indices[2]](2),
132  all_vertices[vertex_indices[3]](2),
133  all_vertices[vertex_indices[4]](2),
134  all_vertices[vertex_indices[5]](2),
135  all_vertices[vertex_indices[6]](2),
136  all_vertices[vertex_indices[7]](2)};
137 
138  /*
139  This is the same Maple script as in the barycenter method above
140  except of that here the shape functions tphi[0]-tphi[7] are ordered
141  according to the lexicographic numbering.
142 
143  x := array(0..7):
144  y := array(0..7):
145  z := array(0..7):
146  tphi[0] := (1-xi)*(1-eta)*(1-zeta):
147  tphi[1] := xi*(1-eta)*(1-zeta):
148  tphi[2] := (1-xi)* eta*(1-zeta):
149  tphi[3] := xi* eta*(1-zeta):
150  tphi[4] := (1-xi)*(1-eta)*zeta:
151  tphi[5] := xi*(1-eta)*zeta:
152  tphi[6] := (1-xi)* eta*zeta:
153  tphi[7] := xi* eta*zeta:
154  x_real := sum(x[s]*tphi[s], s=0..7):
155  y_real := sum(y[s]*tphi[s], s=0..7):
156  z_real := sum(z[s]*tphi[s], s=0..7):
157  with (linalg):
158  J := matrix(3,3, [[diff(x_real, xi), diff(x_real, eta), diff(x_real,
159  zeta)], [diff(y_real, xi), diff(y_real, eta), diff(y_real, zeta)],
160  [diff(z_real, xi), diff(z_real, eta), diff(z_real, zeta)]]):
161  detJ := det (J):
162 
163  measure := simplify ( int ( int ( int (detJ, xi=0..1), eta=0..1),
164  zeta=0..1)):
165 
166  readlib(C):
167 
168  C(measure, optimized);
169 
170  The C code produced by this maple script is further optimized by
171  hand. In particular, division by 12 is performed only once, not
172  hundred of times.
173  */
174 
175  const double t3 = y[3] * x[2];
176  const double t5 = z[1] * x[5];
177  const double t9 = z[3] * x[2];
178  const double t11 = x[1] * y[0];
179  const double t14 = x[4] * y[0];
180  const double t18 = x[5] * y[7];
181  const double t20 = y[1] * x[3];
182  const double t22 = y[5] * x[4];
183  const double t26 = z[7] * x[6];
184  const double t28 = x[0] * y[4];
185  const double t34 =
186  z[3] * x[1] * y[2] + t3 * z[1] - t5 * y[7] + y[7] * x[4] * z[6] +
187  t9 * y[6] - t11 * z[4] - t5 * y[3] - t14 * z[2] + z[1] * x[4] * y[0] -
188  t18 * z[3] + t20 * z[0] - t22 * z[0] - y[0] * x[5] * z[4] - t26 * y[3] +
189  t28 * z[2] - t9 * y[1] - y[1] * x[4] * z[0] - t11 * z[5];
190  const double t37 = y[1] * x[0];
191  const double t44 = x[1] * y[5];
192  const double t46 = z[1] * x[0];
193  const double t49 = x[0] * y[2];
194  const double t52 = y[5] * x[7];
195  const double t54 = x[3] * y[7];
196  const double t56 = x[2] * z[0];
197  const double t58 = x[3] * y[2];
198  const double t64 = -x[6] * y[4] * z[2] - t37 * z[2] + t18 * z[6] -
199  x[3] * y[6] * z[2] + t11 * z[2] + t5 * y[0] +
200  t44 * z[4] - t46 * y[4] - t20 * z[7] - t49 * z[6] -
201  t22 * z[1] + t52 * z[3] - t54 * z[2] - t56 * y[4] -
202  t58 * z[0] + y[1] * x[2] * z[0] + t9 * y[7] + t37 * z[4];
203  const double t66 = x[1] * y[7];
204  const double t68 = y[0] * x[6];
205  const double t70 = x[7] * y[6];
206  const double t73 = z[5] * x[4];
207  const double t76 = x[6] * y[7];
208  const double t90 = x[4] * z[0];
209  const double t92 = x[1] * y[3];
210  const double t95 = -t66 * z[3] - t68 * z[2] - t70 * z[2] + t26 * y[5] -
211  t73 * y[6] - t14 * z[6] + t76 * z[2] - t3 * z[6] +
212  x[6] * y[2] * z[4] - z[3] * x[6] * y[2] + t26 * y[4] -
213  t44 * z[3] - x[1] * y[2] * z[0] + x[5] * y[6] * z[4] +
214  t54 * z[5] + t90 * y[2] - t92 * z[2] + t46 * y[2];
215  const double t102 = x[2] * y[0];
216  const double t107 = y[3] * x[7];
217  const double t114 = x[0] * y[6];
218  const double t125 =
219  y[0] * x[3] * z[2] - z[7] * x[5] * y[6] - x[2] * y[6] * z[4] +
220  t102 * z[6] - t52 * z[6] + x[2] * y[4] * z[6] - t107 * z[5] - t54 * z[6] +
221  t58 * z[6] - x[7] * y[4] * z[6] + t37 * z[5] - t114 * z[4] + t102 * z[4] -
222  z[1] * x[2] * y[0] + t28 * z[6] - y[5] * x[6] * z[4] -
223  z[5] * x[1] * y[4] - t73 * y[7];
224  const double t129 = z[0] * x[6];
225  const double t133 = y[1] * x[7];
226  const double t145 = y[1] * x[5];
227  const double t156 = t90 * y[6] - t129 * y[4] + z[7] * x[2] * y[6] -
228  t133 * z[5] + x[5] * y[3] * z[7] - t26 * y[2] -
229  t70 * z[3] + t46 * y[3] + z[5] * x[7] * y[4] +
230  z[7] * x[3] * y[6] - t49 * z[4] + t145 * z[7] -
231  x[2] * y[7] * z[6] + t70 * z[5] + t66 * z[5] -
232  z[7] * x[4] * y[6] + t18 * z[4] + x[1] * y[4] * z[0];
233  const double t160 = x[5] * y[4];
234  const double t165 = z[1] * x[7];
235  const double t178 = z[1] * x[3];
236  const double t181 =
237  t107 * z[6] + t22 * z[7] + t76 * z[3] + t160 * z[1] - x[4] * y[2] * z[6] +
238  t70 * z[4] + t165 * y[5] + x[7] * y[2] * z[6] - t76 * z[5] - t76 * z[4] +
239  t133 * z[3] - t58 * z[1] + y[5] * x[0] * z[4] + t114 * z[2] - t3 * z[7] +
240  t20 * z[2] + t178 * y[7] + t129 * y[2];
241  const double t207 = t92 * z[7] + t22 * z[6] + z[3] * x[0] * y[2] -
242  x[0] * y[3] * z[2] - z[3] * x[7] * y[2] - t165 * y[3] -
243  t9 * y[0] + t58 * z[7] + y[3] * x[6] * z[2] +
244  t107 * z[2] + t73 * y[0] - x[3] * y[5] * z[7] +
245  t3 * z[0] - t56 * y[6] - z[5] * x[0] * y[4] +
246  t73 * y[1] - t160 * z[6] + t160 * z[0];
247  const double t228 = -t44 * z[7] + z[5] * x[6] * y[4] - t52 * z[4] -
248  t145 * z[4] + t68 * z[4] + t92 * z[5] - t92 * z[0] +
249  t11 * z[3] + t44 * z[0] + t178 * y[5] - t46 * y[5] -
250  t178 * y[0] - t145 * z[0] - t20 * z[5] - t37 * z[3] -
251  t160 * z[7] + t145 * z[3] + x[4] * y[6] * z[2];
252 
253  return (t34 + t64 + t95 + t125 + t156 + t181 + t207 + t228) / 12.;
254  }
255 
256 
257 
258  namespace
259  {
260  // the following class is only
261  // needed in 2d, so avoid trouble
262  // with compilers warning otherwise
263  class Rotate2d
264  {
265  public:
266  explicit Rotate2d(const double angle)
267  : angle(angle)
268  {}
269  Point<2>
270  operator()(const Point<2> &p) const
271  {
272  return {std::cos(angle) * p(0) - std::sin(angle) * p(1),
273  std::sin(angle) * p(0) + std::cos(angle) * p(1)};
274  }
275 
276  private:
277  const double angle;
278  };
279  } // namespace
280 
281 
282 
283  template <>
284  void
286  {
287  transform(Rotate2d(angle), triangulation);
288  }
289 
290 
291 
292  template <>
293  void
295  {
296  (void)angle;
297  (void)triangulation;
298 
299  AssertThrow(
300  false, ExcMessage("GridTools::rotate() is not available for dim = 3."));
301  }
302 } /* namespace GridTools */
303 
std::vector< unsigned int > vertex_indices
Definition: tria.cc:2244
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
double cell_measure< 3 >(const std::vector< Point< 3 >> &all_vertices, const unsigned int(&vertex_indices)[GeometryInfo< 3 >::vertices_per_cell])
const double angle
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
void rotate(const double angle, Triangulation< dim > &triangulation)
double cell_measure< 2 >(const std::vector< Point< 2 >> &all_vertices, const unsigned int(&vertex_indices)[GeometryInfo< 2 >::vertices_per_cell])
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
double cell_measure< 1 >(const std::vector< Point< 1 >> &all_vertices, const unsigned int(&vertex_indices)[GeometryInfo< 1 >::vertices_per_cell])
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation