dune-vtk 2.8
Loading...
Searching...
No Matches
lagrangegridcreator.hh
Go to the documentation of this file.
1#pragma once
2
3#include <cassert>
4#include <cstdint>
5#include <limits>
6#include <optional>
7#include <vector>
8
9#include <dune/common/exceptions.hh>
10#include <dune/common/hybridutilities.hh>
11#include <dune/geometry/utility/typefromvertexcount.hh>
12#include <dune/geometry/multilineargeometry.hh>
13#include <dune/localfunctions/lagrange.hh>
14#include <dune/grid/common/gridfactory.hh>
15
16#include <dune/vtk/types.hh>
20
21namespace Dune
22{
23 namespace Vtk
24 {
25 // \brief Create a grid from data that represents higher (lagrange) cells.
37 template <class GridType>
39 : public GridCreatorInterface<GridType, LagrangeGridCreator<GridType>>
40 {
44
45 using Nodes = std::vector<GlobalCoordinate>;
46
48 {
49 GeometryType type; //< Geometry type of the element
50 std::vector<std::int64_t> nodes; //< Indices of the w.r.t. `nodes_` vector
51 std::vector<unsigned int> corners; //< Insertion-indices of the element corner nodes
52 };
53
54 using Parametrization = std::vector<ElementParametrization>;
55 using Element = typename GridType::template Codim<0>::Entity;
56 using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
57
59 class LocalFunction;
60
61 public:
62 using LocalGeometry = MultiLinearGeometry<typename Element::Geometry::ctype,Element::dimension,Element::dimension>;
63
64 public:
65 using Super::Super;
66 using Super::factory;
67
69 void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
70 std::vector<std::uint64_t> const& /*point_ids*/)
71 {
72 // store point coordinates in member variable
73 nodes_ = points;
74 }
75
76 template <class F>
77 using HasParametrizedElements = decltype(std::declval<F>().insertElement(std::declval<GeometryType>(),
78 std::declval<std::vector<unsigned int> const&>(), std::declval<std::function<GlobalCoordinate(LocalCoordinate)>>()));
79
81 void insertElementsImpl (std::vector<std::uint8_t> const& types,
82 std::vector<std::int64_t> const& offsets,
83 std::vector<std::int64_t> const& connectivity)
84 {
85 assert(nodes_.size() > 0);
86
87 // mapping of node index to element-vertex index
88 std::vector<std::int64_t> elementVertices(nodes_.size(), -1);
89 parametrization_.reserve(types.size());
90
91 std::int64_t vertexIndex = 0;
92 for (std::size_t i = 0; i < types.size(); ++i) {
93 auto type = Vtk::to_geometry(types[i]);
94 if (type.dim() != GridType::dimension)
95 continue;
96
97 Vtk::CellType cellType{type};
98 auto refElem = referenceElement<double,GridType::dimension>(type);
99
100 std::int64_t shift = (i == 0 ? 0 : offsets[i-1]);
101 int nNodes = offsets[i] - shift;
102 int nVertices = refElem.size(GridType::dimension);
103
104 // insert vertices into grid and construct element vertices
105 std::vector<unsigned int> element(nVertices);
106 for (int j = 0; j < nVertices; ++j) {
107 auto index = connectivity.at(shift + j);
108 auto& vertex = elementVertices.at(index);
109 if (vertex < 0) {
110 factory().insertVertex(nodes_.at(index));
111 vertex = vertexIndex++;
112 }
113 element[j] = vertex;
114 }
115
116 // permute element indices
117 if (!cellType.noPermutation()) {
118 // apply index permutation
119 std::vector<unsigned int> cell(element.size());
120 for (std::size_t j = 0; j < element.size(); ++j)
121 cell[j] = element[cellType.permutation(j)];
122 std::swap(element, cell);
123 }
124
125 // fill vector of element parametrizations
126 parametrization_.push_back(ElementParametrization{type});
127 auto& param = parametrization_.back();
128
129 param.nodes.resize(nNodes);
130 for (int j = 0; j < nNodes; ++j)
131 param.nodes[j] = connectivity.at(shift + j);
132 param.corners = element;
133
134 // try to create element with parametrization
135 if constexpr (Std::is_detected_v<HasParametrizedElements, GridFactory<GridType>>) {
136 try {
137 factory().insertElement(type, element,
138 localParametrization(parametrization_.size()-1));
139 } catch (Dune::GridError const& /* notImplemented */) {
140 factory().insertElement(type, element);
141 }
142 } else {
143 factory().insertElement(type, element);
144 }
145 }
146 }
147
149
155 LocalParametrization localParametrization (unsigned int insertionIndex) const
156 {
157 assert(!nodes_.empty() && !parametrization_.empty());
158 auto const& localParam = parametrization_.at(insertionIndex);
159 return LocalParametrization{nodes_, localParam, order(localParam)};
160 }
161
163
169 {
170 VTK_ASSERT(!nodes_.empty() && !parametrization_.empty());
171
172 unsigned int insertionIndex = factory().insertionIndex(element);
173 auto const& localParam = parametrization_.at(insertionIndex);
174 VTK_ASSERT(element.type() == localParam.type);
175
176 return {nodes_, localParam, order(localParam), localGeometry(element, localParam)};
177 }
178
180
187 LocalGeometry localGeometry (Element const& element) const
188 {
189 VTK_ASSERT(!nodes_.empty() && !parametrization_.empty());
190
191 unsigned int insertionIndex = factory().insertionIndex(element);
192 auto const& localParam = parametrization_.at(insertionIndex);
193 VTK_ASSERT(element.type() == localParam.type);
194
195 return localGeometry(element, localParam);
196 }
197
198 private:
199 // implementation details of localGeometry()
200 LocalGeometry localGeometry (Element const& element, ElementParametrization const& localParam) const
201 {
202 // collect indices of vertices
203 std::vector<unsigned int> indices(element.subEntities(GridType::dimension));
204 for (unsigned int i = 0; i < element.subEntities(GridType::dimension); ++i)
205 indices[i] = factory().insertionIndex(element.template subEntity<GridType::dimension>(i));
206
207 // calculate permutation vector
208 std::vector<unsigned int> permutation(indices.size());
209 for (std::size_t i = 0; i < indices.size(); ++i) {
210 auto it = std::find(localParam.corners.begin(), localParam.corners.end(), indices[i]);
211 VTK_ASSERT(it != localParam.corners.end());
212 permutation[i] = std::distance(localParam.corners.begin(), it);
213 }
214
215 auto refElem = referenceElement<typename Element::Geometry::ctype,Element::dimension>(localParam.type);
216 std::vector<LocalCoordinate> corners(permutation.size());
217 for (std::size_t i = 0; i < permutation.size(); ++i)
218 corners[i] = refElem.position(permutation[i], Element::dimension);
219
220 return {localParam.type, corners};
221 }
222
223 public:
224
226 int order (GeometryType type, std::size_t nNodes) const
227 {
228 for (int o = 1; o <= int(nNodes); ++o)
229 if (numLagrangePoints(type, o) == std::size_t(nNodes))
230 return o;
231
232 return 1;
233 }
234
235 int order (ElementParametrization const& localParam) const
236 {
237 return order(localParam.type, localParam.nodes.size());
238 }
239
241 int order () const
242 {
243 assert(!parametrization_.empty());
244 auto const& localParam = parametrization_.front();
245 return order(localParam);
246 }
247
248 public:
250
264 {
265 return LocalFunction{gridCreator};
266 }
267
269 {
270 return LocalFunction{gridCreator};
271 }
272
274 {
275 DUNE_THROW(Dune::Exception, "Cannot pass temporary LagrangeGridCreator to localFunction(). Pass an lvalue-reference instead.");
276 return LocalFunction{gridCreator};
277 }
278
279 std::string name () const
280 {
281 return "LagrangeParametrization";
282 }
283
284 int numComponents () const
285 {
286 return GlobalCoordinate::size();
287 }
288
290 {
291 return dataTypeOf<typename Element::Geometry::ctype>();
292 }
293
295 {
296 using Grid = GridType;
298 };
299
302 {
303 assert(false && "Should not be used!");
304 return EntitySet{};
305 }
306
309 {
310 assert(false && "Should not be used!");
311 return GlobalCoordinate{};
312 }
313
314 private:
316 Nodes nodes_;
317
319 Parametrization parametrization_;
320 };
321
322 // deduction guides
323 template <class Grid>
324 LagrangeGridCreator(GridFactory<Grid>&)
326
327 template <class GridType, class FieldType, class Context>
328 struct AssociatedGridFunction<LagrangeGridCreator<GridType>, FieldType, Context>
329 {
331 };
332
333 template <class Grid>
335 {
336 using ctype = typename Grid::ctype;
337
338 using GlobalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::GlobalCoordinate;
339 using LocalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::LocalCoordinate;
340 using LocalGeometry = MultiLinearGeometry<ctype,Grid::dimension,Grid::dimension>;
341
342 using LocalFE = LagrangeLocalFiniteElement<Vtk::LagrangePointSet, Grid::dimension, ctype, ctype>;
343 using LocalBasis = typename LocalFE::Traits::LocalBasisType;
344 using LocalBasisTraits = typename LocalBasis::Traits;
345
346 public:
348 template <class Nodes, class LocalParam>
349 LocalParametrization (Nodes const& nodes, LocalParam const& param, int order)
350 : localFE_(param.type, order)
351 , localNodes_(param.nodes.size())
352 {
353 for (std::size_t i = 0; i < localNodes_.size(); ++i)
354 localNodes_[i] = nodes[param.nodes[i]];
355 }
356
358 template <class Nodes, class LocalParam, class LG>
359 LocalParametrization (Nodes const& nodes, LocalParam const& param, int order, LG&& localGeometry)
360 : LocalParametrization(nodes, param, order)
361 {
362 localGeometry_.emplace(std::forward<LG>(localGeometry));
363 }
364
366 template <class LocalCoordinate>
367 GlobalCoordinate operator() (LocalCoordinate const& local) const
368 {
369 // map coordinates if element corners are permuted
370 LocalCoordinate x = localGeometry_ ? localGeometry_->global(local) : local;
371
372 LocalBasis const& localBasis = localFE_.localBasis();
373 localBasis.evaluateFunction(x, shapeValues_);
374 assert(shapeValues_.size() == localNodes_.size());
375
376 using field_type = typename LocalBasisTraits::RangeType::field_type;
377
378 GlobalCoordinate out(0);
379 for (std::size_t i = 0; i < shapeValues_.size(); ++i)
380 out.axpy(field_type(shapeValues_[i]), localNodes_[i]);
381
382 return out;
383 }
384
385 private:
386 LocalFE localFE_;
387 std::vector<GlobalCoordinate> localNodes_;
388 std::optional<LocalGeometry> localGeometry_;
389
390 mutable std::vector<typename LocalBasisTraits::RangeType> shapeValues_;
391 };
392
393
394 template <class Grid>
396 {
397 using ctype = typename Grid::ctype;
398 using LocalContext = typename Grid::template Codim<0>::Entity;
399 using GlobalCoordinate = typename LocalContext::Geometry::GlobalCoordinate;
400 using LocalCoordinate = typename LocalContext::Geometry::LocalCoordinate;
401 using LocalParametrization = typename LagrangeGridCreator::LocalParametrization;
402
403 public:
404 explicit LocalFunction (LagrangeGridCreator const& gridCreator)
405 : gridCreator_(&gridCreator)
406 {}
407
408 explicit LocalFunction (LagrangeGridCreator&& gridCreator) = delete;
409
411 void bind (LocalContext const& element)
412 {
413 localContext_ = element;
414 localParametrization_.emplace(gridCreator_->localParametrization(element));
415 }
416
417 void unbind () { /* do nothing */ }
418
420 GlobalCoordinate operator() (LocalCoordinate const& local) const
421 {
422 assert(!!localParametrization_);
423 return (*localParametrization_)(local);
424 }
425
427 LocalContext const& localContext () const
428 {
429 return localContext_;
430 }
431
432 private:
433 LagrangeGridCreator const* gridCreator_;
434
435 LocalContext localContext_;
436 std::optional<LocalParametrization> localParametrization_;
437 };
438
439 } // end namespace Vtk
440} // end namespace Dune
#define VTK_ASSERT(cond)
check if condition cond holds; otherwise, throw a VtkError.
Definition: errors.hh:29
Definition: writer.hh:13
DataTypes
Definition: types.hh:52
GeometryType to_geometry(std::uint8_t cell)
Definition: types.cc:146
Base class for grid creators in a CRTP style.
Definition: gridcreatorinterface.hh:25
GridFactory< Grid > & factory()
Return the associated GridFactory.
Definition: gridcreatorinterface.hh:77
typename Grid::template Codim< 0 >::Entity::Geometry::GlobalCoordinate GlobalCoordinate
Definition: gridcreatorinterface.hh:28
Definition: lagrangegridcreator.hh:40
typename GridType::template Codim< 0 >::Entity Element
Definition: lagrangegridcreator.hh:55
GridFactory< Grid > & factory()
Return the associated GridFactory.
Definition: gridcreatorinterface.hh:77
void insertElementsImpl(std::vector< std::uint8_t > const &types, std::vector< std::int64_t > const &offsets, std::vector< std::int64_t > const &connectivity)
Implementation of the interface function insertElements()
Definition: lagrangegridcreator.hh:81
MultiLinearGeometry< typename Element::Geometry::ctype, Element::dimension, Element::dimension > LocalGeometry
Definition: lagrangegridcreator.hh:62
int numComponents() const
Definition: lagrangegridcreator.hh:284
EntitySet entitySet() const
Dummy function returning a placeholder entityset.
Definition: lagrangegridcreator.hh:301
void insertVerticesImpl(std::vector< GlobalCoordinate > const &points, std::vector< std::uint64_t > const &)
Implementation of the interface function insertVertices()
Definition: lagrangegridcreator.hh:69
decltype(std::declval< F >().insertElement(std::declval< GeometryType >(), std::declval< std::vector< unsigned int > const & >(), std::declval< std::function< GlobalCoordinate(LocalCoordinate)> >())) HasParametrizedElements
Definition: lagrangegridcreator.hh:78
friend LocalFunction localFunction(LagrangeGridCreator &&gridCreator)
Definition: lagrangegridcreator.hh:273
int order(GeometryType type, std::size_t nNodes) const
Determine lagrange order from number of points.
Definition: lagrangegridcreator.hh:226
int order() const
Determine lagrange order from number of points from the first element parametrization.
Definition: lagrangegridcreator.hh:241
LocalGeometry localGeometry(Element const &element) const
Construct a transformation of local element coordinates.
Definition: lagrangegridcreator.hh:187
GlobalCoordinate operator()(GlobalCoordinate const &) const
Dummy function returning a placeholder entityset.
Definition: lagrangegridcreator.hh:308
std::string name() const
Definition: lagrangegridcreator.hh:279
int order(ElementParametrization const &localParam) const
Definition: lagrangegridcreator.hh:235
typename Element::Geometry::LocalCoordinate LocalCoordinate
Definition: lagrangegridcreator.hh:56
LocalParametrization localParametrization(Element const &element) const
Construct an element parametrization.
Definition: lagrangegridcreator.hh:168
LocalParametrization localParametrization(unsigned int insertionIndex) const
Construct an element parametrization.
Definition: lagrangegridcreator.hh:155
std::vector< GlobalCoordinate > Nodes
Definition: lagrangegridcreator.hh:45
friend LocalFunction localFunction(LagrangeGridCreator &gridCreator)
Local function representing the parametrization of the grid.
Definition: lagrangegridcreator.hh:263
Vtk::DataTypes dataType() const
Definition: lagrangegridcreator.hh:289
std::vector< ElementParametrization > Parametrization
Definition: lagrangegridcreator.hh:54
typename Super::GlobalCoordinate GlobalCoordinate
Definition: lagrangegridcreator.hh:43
friend LocalFunction localFunction(LagrangeGridCreator const &gridCreator)
Definition: lagrangegridcreator.hh:268
Definition: lagrangegridcreator.hh:48
std::vector< unsigned int > corners
Definition: lagrangegridcreator.hh:51
GeometryType type
Definition: lagrangegridcreator.hh:49
std::vector< std::int64_t > nodes
Definition: lagrangegridcreator.hh:50
Definition: lagrangegridcreator.hh:295
typename Self::GlobalCoordinate GlobalCoordinate
Definition: lagrangegridcreator.hh:297
GridType Grid
Definition: lagrangegridcreator.hh:296
Definition: lagrangegridcreator.hh:335
LocalParametrization(Nodes const &nodes, LocalParam const &param, int order)
Construct a local element parametrization.
Definition: lagrangegridcreator.hh:349
LocalParametrization(Nodes const &nodes, LocalParam const &param, int order, LG &&localGeometry)
Construct a local element parametrization for elements with permuted corners.
Definition: lagrangegridcreator.hh:359
Definition: lagrangegridcreator.hh:396
LocalFunction(LagrangeGridCreator &&gridCreator)=delete
LocalFunction(LagrangeGridCreator const &gridCreator)
Definition: lagrangegridcreator.hh:404
LocalContext const & localContext() const
Return the bound element.
Definition: lagrangegridcreator.hh:427
void unbind()
Definition: lagrangegridcreator.hh:417
void bind(LocalContext const &element)
Collect a local parametrization on the element.
Definition: lagrangegridcreator.hh:411
Type-Traits to associate a GridFunction to a GridCreator.
Definition: gridfunctions/common.hh:26
Grid-function representing values from a VTK file with local Lagrange interpolation of the values sto...
Definition: lagrangegridfunction.hh:25
Mapping of Dune geometry types to VTK cell types.
Definition: types.hh:160