dune-vtk  0.2
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 
21 namespace 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 DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
230  if (numLagrangePoints(type.id(), type.dim(), o) == std::size_t(nNodes))
231 #else
232  if (numLagrangePoints(type, o) == std::size_t(nNodes))
233 #endif
234  return o;
235 
236  return 1;
237  }
238 
239  int order (ElementParametrization const& localParam) const
240  {
241  return order(localParam.type, localParam.nodes.size());
242  }
243 
245  int order () const
246  {
247  assert(!parametrization_.empty());
248  auto const& localParam = parametrization_.front();
249  return order(localParam);
250  }
251 
252  public:
254 
268  {
269  return LocalFunction{gridCreator};
270  }
271 
272  friend LocalFunction localFunction (LagrangeGridCreator const& gridCreator)
273  {
274  return LocalFunction{gridCreator};
275  }
276 
278  {
279  DUNE_THROW(Dune::Exception, "Cannot pass temporary LagrangeGridCreator to localFunction(). Pass an lvalue-reference instead.");
280  return LocalFunction{gridCreator};
281  }
282 
283  std::string name () const
284  {
285  return "LagrangeParametrization";
286  }
287 
288  int numComponents () const
289  {
290  return GlobalCoordinate::size();
291  }
292 
294  {
295  return dataTypeOf<typename Element::Geometry::ctype>();
296  }
297 
298  struct EntitySet
299  {
300  using Grid = GridType;
302  };
303 
306  {
307  assert(false && "Should not be used!");
308  return EntitySet{};
309  }
310 
313  {
314  assert(false && "Should not be used!");
315  return GlobalCoordinate{};
316  }
317 
318  private:
320  Nodes nodes_;
321 
323  Parametrization parametrization_;
324  };
325 
326  // deduction guides
327  template <class Grid>
328  LagrangeGridCreator(GridFactory<Grid>&)
330 
331  template <class GridType, class FieldType, class Context>
332  struct AssociatedGridFunction<LagrangeGridCreator<GridType>, FieldType, Context>
333  {
335  };
336 
337  template <class Grid>
339  {
340  using ctype = typename Grid::ctype;
341 
342  using GlobalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::GlobalCoordinate;
343  using LocalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::LocalCoordinate;
344  using LocalGeometry = MultiLinearGeometry<ctype,Grid::dimension,Grid::dimension>;
345 
346  using LocalFE = LagrangeLocalFiniteElement<Vtk::LagrangePointSet, Grid::dimension, ctype, ctype>;
347  using LocalBasis = typename LocalFE::Traits::LocalBasisType;
348  using LocalBasisTraits = typename LocalBasis::Traits;
349 
350  public:
352  template <class Nodes, class LocalParam>
353  LocalParametrization (Nodes const& nodes, LocalParam const& param, int order)
354  : localFE_(param.type, order)
355  , localNodes_(param.nodes.size())
356  {
357  for (std::size_t i = 0; i < localNodes_.size(); ++i)
358  localNodes_[i] = nodes[param.nodes[i]];
359  }
360 
362  template <class Nodes, class LocalParam, class LG>
363  LocalParametrization (Nodes const& nodes, LocalParam const& param, int order, LG&& localGeometry)
364  : LocalParametrization(nodes, param, order)
365  {
366  localGeometry_.emplace(std::forward<LG>(localGeometry));
367  }
368 
370  template <class LocalCoordinate>
371  GlobalCoordinate operator() (LocalCoordinate const& local) const
372  {
373  // map coordinates if element corners are permuted
374  LocalCoordinate x = localGeometry_ ? localGeometry_->global(local) : local;
375 
376  LocalBasis const& localBasis = localFE_.localBasis();
377  localBasis.evaluateFunction(x, shapeValues_);
378  assert(shapeValues_.size() == localNodes_.size());
379 
380  using field_type = typename LocalBasisTraits::RangeType::field_type;
381 
382  GlobalCoordinate out(0);
383  for (std::size_t i = 0; i < shapeValues_.size(); ++i)
384  out.axpy(field_type(shapeValues_[i]), localNodes_[i]);
385 
386  return out;
387  }
388 
389  private:
390  LocalFE localFE_;
391  std::vector<GlobalCoordinate> localNodes_;
392  std::optional<LocalGeometry> localGeometry_;
393 
394  mutable std::vector<typename LocalBasisTraits::RangeType> shapeValues_;
395  };
396 
397 
398  template <class Grid>
400  {
401  using ctype = typename Grid::ctype;
402  using LocalContext = typename Grid::template Codim<0>::Entity;
403  using GlobalCoordinate = typename LocalContext::Geometry::GlobalCoordinate;
404  using LocalCoordinate = typename LocalContext::Geometry::LocalCoordinate;
405  using LocalParametrization = typename LagrangeGridCreator::LocalParametrization;
406 
407  public:
408  explicit LocalFunction (LagrangeGridCreator const& gridCreator)
409  : gridCreator_(&gridCreator)
410  {}
411 
412  explicit LocalFunction (LagrangeGridCreator&& gridCreator) = delete;
413 
415  void bind (LocalContext const& element)
416  {
417  localContext_ = element;
418  localParametrization_.emplace(gridCreator_->localParametrization(element));
419  }
420 
421  void unbind () { /* do nothing */ }
422 
424  GlobalCoordinate operator() (LocalCoordinate const& local) const
425  {
426  assert(!!localParametrization_);
427  return (*localParametrization_)(local);
428  }
429 
431  LocalContext const& localContext () const
432  {
433  return localContext_;
434  }
435 
436  private:
437  LagrangeGridCreator const* gridCreator_;
438 
439  LocalContext localContext_;
440  std::optional<LocalParametrization> localParametrization_;
441  };
442 
443  } // end namespace Vtk
444 } // end namespace Dune
#define VTK_ASSERT(cond)
check if condition cond holds; otherwise, throw a VtkError.
Definition: errors.hh:29
Definition: writer.hh:13
LagrangeGridCreator(GridFactory< Grid > &) -> LagrangeGridCreator< Grid >
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
typename Grid::template Codim< 0 >::Entity::Geometry::GlobalCoordinate GlobalCoordinate
Definition: gridcreatorinterface.hh:28
GridFactory< Grid > & factory()
Return the associated GridFactory.
Definition: gridcreatorinterface.hh:77
Definition: lagrangegridcreator.hh:40
typename GridType::template Codim< 0 >::Entity Element
Definition: lagrangegridcreator.hh:55
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:288
EntitySet entitySet() const
Dummy function returning a placeholder entityset.
Definition: lagrangegridcreator.hh:305
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:277
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:245
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:312
std::string name() const
Definition: lagrangegridcreator.hh:283
int order(ElementParametrization const &localParam) const
Definition: lagrangegridcreator.hh:239
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:267
GridFactory< Grid > & factory()
Return the associated GridFactory.
Definition: gridcreatorinterface.hh:77
Vtk::DataTypes dataType() const
Definition: lagrangegridcreator.hh:293
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:272
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:299
typename Self::GlobalCoordinate GlobalCoordinate
Definition: lagrangegridcreator.hh:301
GridType Grid
Definition: lagrangegridcreator.hh:300
Definition: lagrangegridcreator.hh:339
LocalParametrization(Nodes const &nodes, LocalParam const &param, int order)
Construct a local element parametrization.
Definition: lagrangegridcreator.hh:353
LocalParametrization(Nodes const &nodes, LocalParam const &param, int order, LG &&localGeometry)
Construct a local element parametrization for elements with permuted corners.
Definition: lagrangegridcreator.hh:363
Definition: lagrangegridcreator.hh:400
LocalFunction(LagrangeGridCreator &&gridCreator)=delete
LocalFunction(LagrangeGridCreator const &gridCreator)
Definition: lagrangegridcreator.hh:408
void unbind()
Definition: lagrangegridcreator.hh:421
void bind(LocalContext const &element)
Collect a local parametrization on the element.
Definition: lagrangegridcreator.hh:415
LocalContext const & localContext() const
Return the bound element.
Definition: lagrangegridcreator.hh:431
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