Loading...
Searching...
No Matches
Geometry.h
// Copyright (C) 2006-2022 Anders Logg and Garth N. Wells
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier: LGPL-3.0-or-later
#pragma once
#include "Topology.h"
#include <algorithm>
#include <basix/mdspan.hpp>
#include <concepts>
#include <cstdint>
#include <dolfinx/common/IndexMap.h>
#include <dolfinx/common/MPI.h>
#include <dolfinx/common/sort.h>
#include <dolfinx/fem/CoordinateElement.h>
#include <dolfinx/fem/ElementDofLayout.h>
#include <dolfinx/fem/dofmapbuilder.h>
#include <dolfinx/graph/AdjacencyList.h>
#include <dolfinx/graph/partition.h>
#include <functional>
#include <memory>
#include <span>
#include <utility>
#include <vector>
namespace dolfinx::mesh
{
template <std::floating_point T>
class Geometry
{
public:
using value_type = T;
template <typename U, typename V, typename W>
requires std::is_convertible_v<std::remove_cvref_t<U>,
std::vector<std::vector<std::int32_t>>>
and std::is_convertible_v<std::remove_cvref_t<V>,
std::vector<T>>
and std::is_convertible_v<std::remove_cvref_t<W>,
std::vector<std::int64_t>>
Geometry(
std::shared_ptr<const common::IndexMap> index_map, U&& dofmaps,
const std::vector<fem::CoordinateElement<
typename std::remove_reference_t<typename V::value_type>>>& elements,
V&& x, int dim, W&& input_global_indices)
: _dim(dim), _dofmaps(std::forward<U>(dofmaps)),
_index_map(std::move(index_map)), _cmaps(elements),
_x(std::forward<V>(x)),
_input_global_indices(std::forward<W>(input_global_indices))
{
assert(_x.size() % 3 == 0);
if (_x.size() / 3 != _input_global_indices.size())
throw std::runtime_error("Geometry size mismatch.");
if (_dofmaps.size() != _cmaps.size())
{
throw std::runtime_error("Geometry number of dofmaps not equal to the "
"number of coordinate elements.");
}
// TODO: check that elements dim == number of dofmap columns
}
Geometry(const Geometry&) = default;
Geometry(Geometry&&) = default;
~Geometry() = default;
Geometry& operator=(const Geometry&) = delete;
Geometry& operator=(Geometry&&) = default;
int dim() const { return _dim; }
md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>> dofmap() const
{
if (_dofmaps.size() != 1)
throw std::runtime_error("Multiple dofmaps");
return this->dofmap(0);
}
md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>>
dofmap(std::size_t i) const
{
std::size_t ndofs = _cmaps.at(i).dim();
return md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>>(
_dofmaps.at(i).data(), _dofmaps.at(i).size() / ndofs, ndofs);
}
std::shared_ptr<const common::IndexMap> index_map() const
{
return _index_map;
}
std::span<const value_type> x() const { return _x; }
std::span<value_type> x() { return _x; }
const std::vector<fem::CoordinateElement<value_type>>& cmaps() const
{
return _cmaps;
}
const fem::CoordinateElement<value_type>& cmap() const
{
if (_cmaps.size() > 1)
throw std::runtime_error("Multiple cmaps.");
return _cmaps.front();
}
const std::vector<std::int64_t>& input_global_indices() const
{
return _input_global_indices;
}
private:
// Geometric dimension
int _dim;
// Map per cell for extracting coordinate data for each cmap
std::vector<std::vector<std::int32_t>> _dofmaps;
// IndexMap for geometry 'dofmap'
std::shared_ptr<const common::IndexMap> _index_map;
// The coordinate elements
std::vector<fem::CoordinateElement<value_type>> _cmaps;
// Coordinates for all points stored as a contiguous array (row-major,
// column size = 3)
std::vector<value_type> _x;
// Global indices as provided on Geometry creation
std::vector<std::int64_t> _input_global_indices;
};
template <typename U, typename V, typename W>
Geometry(std::shared_ptr<const common::IndexMap>, U&&,
const std::vector<fem::CoordinateElement<
typename std::remove_reference_t<typename V::value_type>>>&,
V&&, int, W&&)
-> Geometry<typename std::remove_cvref_t<typename V::value_type>>;
template <typename U>
Geometry<typename std::remove_reference_t<typename U::value_type>>
create_geometry(const Topology& topology,
const std::vector<fem::CoordinateElement<
std::remove_reference_t<typename U::value_type>>>& elements,
std::span<const std::int64_t> nodes,
std::span<const std::int64_t> xdofs, const U& x, int dim,
const std::function<std::vector<int>(
const graph::AdjacencyList<std::int32_t>&)>& reorder_fn
= nullptr)
{
spdlog::info("Create Geometry (multiple)");
assert(std::ranges::is_sorted(nodes));
using T = typename std::remove_reference_t<typename U::value_type>;
// Check elements match cell types in topology
const int tdim = topology.dim();
const std::size_t num_cell_types = topology.entity_types(tdim).size();
if (elements.size() != num_cell_types)
throw std::runtime_error("Mismatch between topology and geometry.");
std::vector<fem::ElementDofLayout> dof_layouts;
dof_layouts.reserve(elements.size());
for (auto& el : elements)
dof_layouts.push_back(el.create_dof_layout());
spdlog::info("Got {} dof layouts", dof_layouts.size());
// Build 'geometry' dofmap on the topology
auto [_dof_index_map, bs, dofmaps]
= fem::build_dofmap_data(topology.index_maps(topology.dim())[0]->comm(),
topology, dof_layouts, reorder_fn);
auto dof_index_map
= std::make_shared<common::IndexMap>(std::move(_dof_index_map));
// If the mesh has higher order geometry, permute the dofmap
if (elements.front().needs_dof_permutations())
{
const std::int32_t num_cells
= topology.connectivity(topology.dim(), 0)->num_nodes();
const std::vector<std::uint32_t>& cell_info
= topology.get_cell_permutation_info();
int d = elements.front().dim();
for (std::int32_t cell = 0; cell < num_cells; ++cell)
{
std::span dofs(dofmaps.front().data() + cell * d, d);
elements.front().permute_inv(dofs, cell_info[cell]);
}
}
spdlog::info("Calling compute_local_to_global");
// Compute local-to-global map from local indices in dofmap to the
// corresponding global indices in cells, and pass to function to
// compute local (dof) to local (position in coords) map from (i)
// local-to-global for dofs and (ii) local-to-global for entries in
// coords
spdlog::info("xdofs.size = {}", xdofs.size());
std::vector<std::int32_t> all_dofmaps;
std::stringstream s;
for (auto q : dofmaps)
{
s << q.size() << " ";
all_dofmaps.insert(all_dofmaps.end(), q.begin(), q.end());
}
spdlog::info("dofmap sizes = {}", s.str());
spdlog::info("all_dofmaps.size = {}", all_dofmaps.size());
spdlog::info("nodes.size = {}", nodes.size());
const std::vector<std::int32_t> l2l = graph::build::compute_local_to_local(
graph::build::compute_local_to_global(xdofs, all_dofmaps), nodes);
// Allocate space for input global indices and copy data
std::vector<std::int64_t> igi(nodes.size());
std::ranges::transform(l2l, igi.begin(),
[&nodes](auto index) { return nodes[index]; });
// Build coordinate dof array, copying coordinates to correct position
assert(x.size() % dim == 0);
const std::size_t shape0 = x.size() / dim;
const std::size_t shape1 = dim;
std::vector<T> xg(3 * shape0, 0);
for (std::size_t i = 0; i < shape0; ++i)
{
std::copy_n(std::next(x.begin(), shape1 * l2l[i]), shape1,
std::next(xg.begin(), 3 * i));
}
spdlog::info("Creating geometry with {} dofmaps", dof_layouts.size());
return Geometry(dof_index_map, std::move(dofmaps), elements, std::move(xg),
dim, std::move(igi));
}
} // namespace dolfinx::mesh
dolfinx::fem::CoordinateElement
Definition CoordinateElement.h:38
dolfinx::graph::AdjacencyList
This class provides a static adjacency list data structure.
Definition AdjacencyList.h:38
dolfinx::mesh::Geometry
Geometry stores the geometry imposed on a mesh.
Definition Geometry.h:34
dolfinx::mesh::Geometry::~Geometry
~Geometry()=default
Destructor.
dolfinx::mesh::Geometry::x
std::span< value_type > x()
Access geometry degrees-of-freedom data (non-const version).
Definition Geometry.h:151
dolfinx::mesh::Geometry::Geometry
Geometry(std::shared_ptr< const common::IndexMap > index_map, U &&dofmaps, const std::vector< fem::CoordinateElement< typename std::remove_reference_t< typename V::value_type > > > &elements, V &&x, int dim, W &&input_global_indices)
Constructor of object that holds mesh geometry data.
Definition Geometry.h:69
dolfinx::mesh::Geometry::dofmap
md::mdspan< const std::int32_t, md::dextents< std::size_t, 2 > > dofmap(std::size_t i) const
Degree-of-freedom map associated with the ith coordinate map element in the geometry.
Definition Geometry.h:126
dolfinx::mesh::Geometry::index_map
std::shared_ptr< const common::IndexMap > index_map() const
Index map for the geometry 'degrees-of-freedom'.
Definition Geometry.h:135
dolfinx::mesh::Geometry::Geometry
Geometry(Geometry &&)=default
Move constructor.
dolfinx::mesh::Geometry::Geometry
Geometry(const Geometry &)=default
Copy constructor.
dolfinx::mesh::Geometry::cmap
const fem::CoordinateElement< value_type > & cmap() const
The element that describes the geometry map.
Definition Geometry.h:167
dolfinx::mesh::Geometry::dim
int dim() const
Return dimension of the Euclidean coordinate system.
Definition Geometry.h:108
dolfinx::mesh::Geometry::operator=
Geometry & operator=(Geometry &&)=default
Move assignment.
dolfinx::mesh::Geometry::dofmap
md::mdspan< const std::int32_t, md::dextents< std::size_t, 2 > > dofmap() const
DofMap for the geometry.
Definition Geometry.h:112
dolfinx::mesh::Geometry::input_global_indices
const std::vector< std::int64_t > & input_global_indices() const
Global user indices.
Definition Geometry.h:175
dolfinx::mesh::Geometry::x
std::span< const value_type > x() const
Access geometry degrees-of-freedom data (const version).
Definition Geometry.h:144
dolfinx::mesh::Geometry::operator=
Geometry & operator=(const Geometry &)=delete
Copy assignment.
dolfinx::mesh::Geometry::cmaps
const std::vector< fem::CoordinateElement< value_type > > & cmaps() const
The elements that describes the geometry map.
Definition Geometry.h:159
dolfinx::mesh::Geometry::value_type
T value_type
Value type.
Definition Geometry.h:37
dolfinx::mesh::Topology
Topology stores the topology of a mesh, consisting of mesh entities and connectivity (incidence relat...
Definition Topology.h:41
dolfinx::mesh::Topology::get_cell_permutation_info
const std::vector< std::uint32_t > & get_cell_permutation_info() const
Returns the permutation information.
Definition Topology.cpp:860
dolfinx::mesh::Topology::connectivity
std::shared_ptr< const graph::AdjacencyList< std::int32_t > > connectivity(std::array< int, 2 > d0, std::array< int, 2 > d1) const
Get the connectivity from entities of topological dimension d0 to dimension d1.
Definition Topology.cpp:840
dolfinx::mesh::Topology::entity_types
const std::vector< CellType > & entity_types(int dim) const
Entity types in the topology for a given dimension.
Definition Topology.cpp:791
dolfinx::mesh::Topology::dim
int dim() const noexcept
Topological dimension of the mesh.
Definition Topology.cpp:786
dolfinx::mesh::Topology::index_maps
std::vector< std::shared_ptr< const common::IndexMap > > index_maps(int dim) const
Get the index maps that described the parallel distribution of the mesh entities of a given topologic...
Definition Topology.cpp:812
dolfinx::fem::build_dofmap_data
std::tuple< common::IndexMap, int, std::vector< std::vector< std::int32_t > > > build_dofmap_data(MPI_Comm comm, const mesh::Topology &topology, const std::vector< ElementDofLayout > &element_dof_layouts, const std::function< std::vector< int >(const graph::AdjacencyList< std::int32_t > &)> &reorder_fn)
Definition dofmapbuilder.cpp:627
dolfinx::graph::build::compute_local_to_global
std::vector< std::int64_t > compute_local_to_global(std::span< const std::int64_t > global, std::span< const std::int32_t > local)
Definition partition.cpp:534
dolfinx::graph::build::compute_local_to_local
std::vector< std::int32_t > compute_local_to_local(std::span< const std::int64_t > local0_to_global, std::span< const std::int64_t > local1_to_global)
Compute a local0-to-local1 map from two local-to-global maps with common global indices.
Definition partition.cpp:556
dolfinx::mesh
Mesh data structures and algorithms on meshes.
Definition DofMap.h:32
dolfinx::mesh::create_geometry
Geometry< typename std::remove_reference_t< typename U::value_type > > create_geometry(const Topology &topology, const std::vector< fem::CoordinateElement< std::remove_reference_t< typename U::value_type > > > &elements, std::span< const std::int64_t > nodes, std::span< const std::int64_t > xdofs, const U &x, int dim, const std::function< std::vector< int >(const graph::AdjacencyList< std::int32_t > &)> &reorder_fn=nullptr)
Build Geometry from input data.
Definition Geometry.h:237