-
Notifications
You must be signed in to change notification settings - Fork 707
/
mapping_manifold.h
460 lines (383 loc) · 15.3 KB
/
mapping_manifold.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 2021 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------
#ifndef dealii_mapping_manifold_h
#define dealii_mapping_manifold_h
#include <deal.II/base/config.h>
#include <deal.II/base/array_view.h>
#include <deal.II/base/derivative_form.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/fe/mapping.h>
#include <cmath>
DEAL_II_NAMESPACE_OPEN
template <int, int>
class MappingQ;
/*!@addtogroup mapping */
/*@{*/
/**
* This class implements the functionality for Manifold conforming
* mappings. This Mapping computes the transformation between the
* reference and real cell by exploiting the geometrical information
* coming from the underlying Manifold object.
*
* Quadrature points computed using this mapping lie on the exact
* geometrical objects, and tangent and normal vectors computed using
* this class are tangent and normal to the underlying geometry. This
* is in contrast with the MappingQ class, which approximates the
* geometry using a polynomial of some order, and then computes the
* normals and tangents using the approximated surface.
*
* @warning It is not possible, for mathematical reasons, for one to use this
* class with a geometry described by a SphericalManifold: see the note in
* that class for more information.
*/
template <int dim, int spacedim = dim>
class MappingManifold : public Mapping<dim, spacedim>
{
public:
/**
* Constructor.
*/
MappingManifold() = default;
/**
* Copy constructor.
*/
MappingManifold(const MappingManifold<dim, spacedim> &mapping);
// for documentation, see the Mapping base class
virtual std::unique_ptr<Mapping<dim, spacedim>>
clone() const override;
/**
* Always returns @p true because this class assumes that the
* vertices always lies on the underlying Manifold.
*/
virtual bool
preserves_vertex_locations() const override;
virtual bool
is_compatible_with(const ReferenceCell &cell_type) const override;
/**
* @name Mapping points between reference and real cells
* @{
*/
// for documentation, see the Mapping base class
virtual Point<spacedim>
transform_unit_to_real_cell(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const Point<dim> &p) const override;
// for documentation, see the Mapping base class
virtual Point<dim>
transform_real_to_unit_cell(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const Point<spacedim> &p) const override;
/**
* @}
*/
/**
* @name Functions to transform tensors from reference to real coordinates
* @{
*/
// for documentation, see the Mapping base class
virtual void
transform(const ArrayView<const Tensor<1, dim>> & input,
const MappingKind kind,
const typename Mapping<dim, spacedim>::InternalDataBase &internal,
const ArrayView<Tensor<1, spacedim>> &output) const override;
// for documentation, see the Mapping base class
virtual void
transform(const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
const MappingKind kind,
const typename Mapping<dim, spacedim>::InternalDataBase &internal,
const ArrayView<Tensor<2, spacedim>> &output) const override;
// for documentation, see the Mapping base class
virtual void
transform(const ArrayView<const Tensor<2, dim>> & input,
const MappingKind kind,
const typename Mapping<dim, spacedim>::InternalDataBase &internal,
const ArrayView<Tensor<2, spacedim>> &output) const override;
// for documentation, see the Mapping base class
virtual void
transform(const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
const MappingKind kind,
const typename Mapping<dim, spacedim>::InternalDataBase &internal,
const ArrayView<Tensor<3, spacedim>> &output) const override;
// for documentation, see the Mapping base class
virtual void
transform(const ArrayView<const Tensor<3, dim>> & input,
const MappingKind kind,
const typename Mapping<dim, spacedim>::InternalDataBase &internal,
const ArrayView<Tensor<3, spacedim>> &output) const override;
/**
* @}
*/
/**
* @name Interface with FEValues
* @{
*/
/**
* Storage for internal data of polynomial mappings. See
* Mapping::InternalDataBase for an extensive description.
*
* For the current class, the InternalData class stores data that is
* computed once when the object is created (in get_data()) as well as data
* the class wants to store from between the call to fill_fe_values(),
* fill_fe_face_values(), or fill_fe_subface_values() until possible later
* calls from the finite element to functions such as transform(). The
* latter class of member variables are marked as 'mutable'.
*/
class InternalData : public Mapping<dim, spacedim>::InternalDataBase
{
public:
/**
* Constructor.
*/
InternalData() = default;
/**
* Initialize the object's member variables related to cell data based on
* the given arguments.
*
* The function also calls compute_shape_function_values() to actually set
* the member variables related to the values and derivatives of the
* mapping shape functions.
*/
void
initialize(const UpdateFlags update_flags,
const Quadrature<dim> &quadrature,
const unsigned int n_original_q_points);
/**
* Initialize the object's member variables related to cell and face data
* based on the given arguments. In order to initialize cell data, this
* function calls initialize().
*/
void
initialize_face(const UpdateFlags update_flags,
const Quadrature<dim> &quadrature,
const unsigned int n_original_q_points);
/**
* Compute the weights associated to the Manifold object, that
* need to be passed when computing the location of the quadrature
* points.
*/
void
compute_manifold_quadrature_weights(const Quadrature<dim> &quadrature);
/**
* Store vertices internally.
*/
void
store_vertices(
const typename Triangulation<dim, spacedim>::cell_iterator &cell) const;
/**
* Return an estimate (in bytes) for the memory consumption of this object.
*/
virtual std::size_t
memory_consumption() const override;
/**
* The current cell vertices.
*
* Computed each.
*/
mutable std::vector<Point<spacedim>> vertices;
/**
* The current cell.
*
* Computed each.
*/
mutable typename Triangulation<dim, spacedim>::cell_iterator cell;
/**
* The actual quadrature on the reference cell.
*
* Computed once.
*/
Quadrature<dim> quad;
/**
* Values of quadrature weights for manifold quadrature
* formulas.
*
* The Manifold class has a function (Manifold::get_new_point())
* that returns new points according to a weighted average of some
* surrounding points on the Manifold. For each quadrature point,
* we call this function with a Quadrature formula constructed
* using the vertices of the current cell, and the values of the
* basis functions of an FE_Q(1) finite element evaluated at the
* quadrature point itself. While the vertices of the cell change
* for every cell, the weights can be computed once for each
* quadrature point. We store this information in the following
* variable, where the first index runs through the quadrature
* points, and the second index runs through the vertex indices.
*
* Computed once.
*/
std::vector<std::vector<double>> cell_manifold_quadrature_weights;
/**
* A vector of weights for use in Manifold::get_new_point(). For
* each point (interior to a cell), we compute the weight each
* vertex has for this point. If the point lies at a vertex, then
* this vertex has weight one and all others have weight zero. If
* the point lies interior to a cell, then the weight every vertex
* has is just the $d$-linear shape functions associated with each
* vertex evaluated at that point.
*
* This array has size GeometryInfo<dim>::vertices_per_cell, but it
* can't be converted into a fixed size array because it is used
* as input for Manifold::get_new_point() which wants to see a
* std::vector<double> for the weights.
*/
mutable std::vector<double> vertex_weights;
/**
* Unit tangential vectors. Used for the computation of boundary forms and
* normal vectors.
*
* This array has `(dim-1) * GeometryInfo::faces_per_cell` entries. The
* first GeometryInfo::faces_per_cell contain the vectors in the first
* tangential direction for each face; the second set of
* GeometryInfo<dim>::faces_per_cell entries contain the vectors in the
* second tangential direction (only in 3d, since there we have 2 tangential
* directions per face), etc.
*
* Filled once.
*/
std::array<std::vector<Tensor<1, dim>>,
GeometryInfo<dim>::faces_per_cell *(dim - 1)>
unit_tangentials;
/**
* Tensors of covariant transformation at each of the quadrature points.
* The matrix stored is the Jacobian * G^{-1}, where G = Jacobian^{t} *
* Jacobian, is the first fundamental form of the map; if dim=spacedim
* then it reduces to the transpose of the inverse of the Jacobian matrix,
* which itself is stored in the @p contravariant field of this structure.
*
* Computed on each cell.
*/
mutable std::vector<DerivativeForm<1, dim, spacedim>> covariant;
/**
* Tensors of contravariant transformation at each of the quadrature
* points. The contravariant matrix is the Jacobian of the transformation,
* i.e. $J_{ij}=dx_i/d\hat x_j$.
*
* Computed on each cell.
*/
mutable std::vector<DerivativeForm<1, dim, spacedim>> contravariant;
/**
* Auxiliary vectors for internal use.
*/
mutable std::vector<std::vector<Tensor<1, spacedim>>> aux;
/**
* The determinant of the Jacobian in each quadrature point. Filled if
* #update_volume_elements.
*/
mutable std::vector<double> volume_elements;
/**
* A pointer to the Manifold in use.
*
* Updated each.
*/
mutable SmartPointer<const Manifold<dim, spacedim>> manifold;
};
private:
// documentation can be found in Mapping::requires_update_flags()
virtual UpdateFlags
requires_update_flags(const UpdateFlags update_flags) const override;
// documentation can be found in Mapping::get_data()
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
get_data(const UpdateFlags, const Quadrature<dim> &quadrature) const override;
using Mapping<dim, spacedim>::get_face_data;
// documentation can be found in Mapping::get_face_data()
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
get_face_data(const UpdateFlags flags,
const hp::QCollection<dim - 1> &quadrature) const override;
// documentation can be found in Mapping::get_subface_data()
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
get_subface_data(const UpdateFlags flags,
const Quadrature<dim - 1> &quadrature) const override;
// documentation can be found in Mapping::fill_fe_values()
virtual CellSimilarity::Similarity
fill_fe_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const CellSimilarity::Similarity cell_similarity,
const Quadrature<dim> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
using Mapping<dim, spacedim>::fill_fe_face_values;
// documentation can be found in Mapping::fill_fe_face_values()
virtual void
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
const hp::QCollection<dim - 1> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
// documentation can be found in Mapping::fill_fe_subface_values()
virtual void
fill_fe_subface_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int subface_no,
const Quadrature<dim - 1> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
/**
* @}
*/
};
/*@}*/
/*----------------------------------------------------------------------*/
#ifndef DOXYGEN
template <int dim, int spacedim>
inline void
MappingManifold<dim, spacedim>::InternalData::store_vertices(
const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
{
vertices.resize(GeometryInfo<dim>::vertices_per_cell);
for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
vertices[i] = cell->vertex(i);
this->cell = cell;
}
template <int dim, int spacedim>
inline void
MappingManifold<dim, spacedim>::InternalData::
compute_manifold_quadrature_weights(const Quadrature<dim> &quad)
{
cell_manifold_quadrature_weights.resize(
quad.size(), std::vector<double>(GeometryInfo<dim>::vertices_per_cell));
for (unsigned int q = 0; q < quad.size(); ++q)
{
for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
{
cell_manifold_quadrature_weights[q][i] =
GeometryInfo<dim>::d_linear_shape_function(quad.point(q), i);
}
}
}
template <int dim, int spacedim>
inline bool
MappingManifold<dim, spacedim>::preserves_vertex_locations() const
{
return true;
}
template <int dim, int spacedim>
bool
MappingManifold<dim, spacedim>::is_compatible_with(
const ReferenceCell &cell_type) const
{
if (cell_type.get_dimension() != dim)
return false; // TODO: or is this an error?
if (cell_type.is_hyper_cube())
return true;
return false;
}
#endif // DOXYGEN
/* -------------- declaration of explicit specializations ------------- */
DEAL_II_NAMESPACE_CLOSE
#endif