preCICE v3.1.1
Loading...
Searching...
No Matches
RadialGeoMultiscaleMapping.cpp
Go to the documentation of this file.
2#include <Eigen/Core>
3#include <algorithm>
4#include <numeric>
5#include "mesh/Mesh.hpp"
6
7namespace precice::mapping {
8
20
22{
23 PRECICE_TRACE(output()->nVertices());
24
25 PRECICE_ASSERT(input().get() != nullptr);
26 PRECICE_ASSERT(output().get() != nullptr);
27
28 size_t const inSize = input()->nVertices();
29 size_t const outSize = output()->nVertices();
30
31 int effectiveCoordinate = static_cast<std::underlying_type_t<MultiscaleType>>(_axis);
32 PRECICE_ASSERT(effectiveCoordinate == static_cast<std::underlying_type_t<MultiscaleType>>(MultiscaleAxis::X) ||
33 effectiveCoordinate == static_cast<std::underlying_type_t<MultiscaleType>>(MultiscaleAxis::Y) ||
34 effectiveCoordinate == static_cast<std::underlying_type_t<MultiscaleType>>(MultiscaleAxis::Z),
35 "Unknown multiscale axis type.");
36
37 if (getConstraint() == CONSISTENT) {
38 PRECICE_DEBUG("Compute consistent mapping");
40 // Determine principle axis midpoints as borders to assign 3D vertices to their respective 1D vertex
41 Eigen::VectorXd axisMidpoints(inSize);
42 auto & inputVerticesRef = input()->vertices();
43 // Order the vertices of the 1D input mesh, to correctly compute the midpoints of the segments between neighboring vertices
44 std::vector<size_t> ordered_vertex_indices(input()->nVertices());
45 std::iota(ordered_vertex_indices.begin(), ordered_vertex_indices.end(), 0);
46 std::stable_sort(ordered_vertex_indices.begin(), ordered_vertex_indices.end(),
47 [effectiveCoordinate, &inputVerticesRef](const size_t aindex, const size_t bindex) {
48 return inputVerticesRef[aindex].coord(effectiveCoordinate) < inputVerticesRef[bindex].coord(effectiveCoordinate);
49 });
50 // Compute the midpoints of the 1D input mesh
51 for (size_t i = 0; i < (inSize - 1); i++) {
52 auto axisPositionCurrent = inputVerticesRef[ordered_vertex_indices[i]].coord(effectiveCoordinate);
53 auto axisPositionNext = inputVerticesRef[ordered_vertex_indices[i] + 1].coord(effectiveCoordinate);
54 axisMidpoints(i) = (axisPositionCurrent + axisPositionNext) / 2;
55 }
56 axisMidpoints(inSize - 1) = std::numeric_limits<double>::max(); // large number, such that vertices after the last midpoint are still assigned
57
58 /*
59 3D vertices are projected onto the 1D axis and the data is then mapped
60 to the nearest neighbors of the 1D vertices in projection space.
61 */
63 _vertexIndicesSpread.reserve(output()->nVertices());
64 for (size_t i = 0; i < outSize; i++) {
65 auto vertexCoord = output()->vertex(i).coord(effectiveCoordinate);
66 size_t index = 0;
67 while (vertexCoord > axisMidpoints(index)) {
68 PRECICE_ASSERT(index + 1 < inSize);
69 ++index;
70 }
72 }
73 } else {
75 /*
76 3D vertices are projected onto the 1D axis and the data is then mapped
77 to (and averaged at) the nearest 1D vertex in projection space.
78 */
79
80 // determine principle axis midpoints as borders to assign 3D vertices to their respective 1D vertex
81 Eigen::VectorXd axisMidpoints(outSize);
82 auto & outputVerticesRef = output()->vertices();
83
84 // Order the vertices of the 1D output mesh, to correctly compute the midpoints of the segments between neighboring vertices
85 std::vector<size_t> ordered_vertex_indices(output()->nVertices());
86 std::iota(ordered_vertex_indices.begin(), ordered_vertex_indices.end(), 0);
87 std::stable_sort(ordered_vertex_indices.begin(), ordered_vertex_indices.end(),
88 [effectiveCoordinate, &outputVerticesRef](const size_t aindex, const size_t bindex) {
89 return outputVerticesRef[aindex].coord(effectiveCoordinate) < outputVerticesRef[bindex].coord(effectiveCoordinate);
90 });
91 // Compute the midpoints of the 1D output mesh
92 for (size_t i = 0; i < (outSize - 1); i++) {
93 auto axisPositionCurrent = output()->vertex(i).coord(effectiveCoordinate);
94 auto axisPositionNext = output()->vertex(i + 1).coord(effectiveCoordinate);
95 axisMidpoints(i) = (axisPositionCurrent + axisPositionNext) / 2;
96 }
97 axisMidpoints(outSize - 1) = std::numeric_limits<double>::max(); // large number, such that vertices after the last midpoint are still assigned
98
99 std::vector<size_t> counters(outSize); // counts number of vertices in between midpoints for averaging
100
101 // Identify which vertex (index) of the 3D mesh corresponds to which vertex (index) of the 1D meesh
103 _vertexIndicesCollect.reserve(input()->nVertices());
104 for (size_t i = 0; i < inSize; i++) {
105 auto vertexCoords = input()->vertex(i).coord(effectiveCoordinate);
106 size_t index = 0;
107 while (vertexCoords > axisMidpoints(index)) {
108 PRECICE_ASSERT(index + 1 < outSize);
109 ++index;
110 }
112 counters[index] += 1;
113 }
114 _vertexCounter = std::move(counters);
115 }
116
117 } else {
119 PRECICE_ASSERT(false, "Not yet implemented");
120 PRECICE_DEBUG("Compute conservative mapping");
121 }
122 _hasComputedMapping = true;
123}
124
132
133void RadialGeoMultiscaleMapping::mapConservative(const time::Sample &inData, Eigen::VectorXd &outData)
134{
136 PRECICE_ASSERT(false, "Not yet implemented");
137 PRECICE_DEBUG("Map conservative");
138}
139
140void RadialGeoMultiscaleMapping::mapConsistent(const time::Sample &inData, Eigen::VectorXd &outData)
141{
143
144 const int inDataDimensions = inData.dataDims;
145 const Eigen::VectorXd &inputValues = inData.values;
146 Eigen::VectorXd & outputValues = outData;
147
148 size_t const inSize = input()->nVertices();
149 size_t const outSize = output()->nVertices();
150
151 PRECICE_ASSERT(!output()->empty());
152 auto outDataDimensions = outputValues.size() / outSize;
153
154 // Check that the number of values for the input and output is right according to their dimensions
155 PRECICE_ASSERT((inputValues.size() / static_cast<std::size_t>(inDataDimensions) == input()->nVertices()),
156 inputValues.size(), inDataDimensions, input()->nVertices());
157 PRECICE_ASSERT((outputValues.size() / outDataDimensions == output()->nVertices()),
158 outputValues.size(), outDataDimensions, output()->nVertices());
159
160 // We currently don't support 1D data, so we need that the user specifies data of the same dimensions on both sides
161 PRECICE_ASSERT(static_cast<std::size_t>(inDataDimensions) == outDataDimensions);
162
163 PRECICE_DEBUG("Map consistent");
165 // assign 1D vertex value to all 3D vertices in vicinity
166 for (size_t i = 0; i < outSize; i++) {
167 outputValues(i * outDataDimensions) = inputValues(_vertexIndicesSpread[i] * inDataDimensions);
168 }
169 } else {
171 /*
172 3D vertices are projected onto the 1D axis and the data is then mapped
173 to (and averaged at) the nearest 1D vertex in projection space.
174 */
175
176 for (size_t i = 0; i < outSize; i++) {
177 outputValues(i * outDataDimensions) = 0;
178 }
179 // assign the 1D vertex the average of all 3D vertex values in vicinity
180 for (size_t i = 0; i < inSize; i++) {
181 outputValues(_vertexIndicesCollect[i] * outDataDimensions) += inputValues(i * inDataDimensions);
182 }
183 for (size_t i = 0; i < outSize; i++) {
184 outputValues(i * outDataDimensions) = outputValues(i * outDataDimensions) / _vertexCounter[i];
185 }
186 }
187}
188
190{
192
194
195 if (getConstraint() == CONSISTENT) {
196 PRECICE_ASSERT(_type == MultiscaleType::SPREAD, "Not yet implemented");
197
198 // tag all vertices of the 1D mesh
199 input()->tagAll();
200
201 } else {
203 PRECICE_ASSERT(false, "Not yet implemented");
204 }
205
206 clear();
207}
208
210{
212 // no operation needed here for the moment
213}
214
216{
217 return "radial-geomultiscale";
218}
219
220} // namespace precice::mapping
#define PRECICE_DEBUG(...)
Definition LogMacros.hpp:64
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
unsigned int index
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
T begin(T... args)
Abstract base class for mapping of data from one mesh to another.
Definition Mapping.hpp:15
mesh::PtrMesh output() const
Returns pointer to output mesh.
Definition Mapping.cpp:91
Constraint
Specifies additional constraints for a mapping.
Definition Mapping.hpp:29
mesh::PtrMesh input() const
Returns pointer to input mesh.
Definition Mapping.cpp:86
bool _hasComputedMapping
Flag to indicate whether computeMapping() has been called.
Definition Mapping.hpp:208
Constraint getConstraint() const
Returns the constraint (consistent/conservative) of the mapping.
Definition Mapping.cpp:45
void setInputRequirement(MeshRequirement requirement)
Sets the mesh requirement for the input mesh.
Definition Mapping.cpp:96
void setOutputRequirement(MeshRequirement requirement)
Sets the mesh requirement for the output mesh.
Definition Mapping.cpp:102
InitialGuessRequirement
Specifies whether the mapping requires an initial guess.
Definition Mapping.hpp:63
MultiscaleType
Geometric multiscale type of the mapping (spread or collect).
RadialGeoMultiscaleMapping(Constraint constraint, int dimensions, MultiscaleType type, MultiscaleAxis axis)
Constructor.
std::vector< size_t > _vertexCounter
counts number of vertices between midpoints for averaging
MultiscaleAxis _axis
main axis along which radial geometric multiscale coupling happens
void tagMeshFirstRound() override
Method used by partition. Tags vertices that could be owned by this rank.
std::vector< size_t > _vertexIndicesSpread
computed vertex indices to map data from input vertices to output vertices and vice versa
std::string getName() const final override
Returns name of the mapping.
MultiscaleType _type
type of mapping, namely spread or collect
void clear() override
Removes a computed mapping.
void tagMeshSecondRound() override
Method used by partition. Tags vertices that can be filtered out.
void mapConservative(const time::Sample &inData, Eigen::VectorXd &outData) override
Maps data using a conservative constraint.
void computeMapping() override
Takes care of compute-heavy operations needed only once to set up the mapping.
void mapConsistent(const time::Sample &inData, Eigen::VectorXd &outData) override
Maps data using a consistent constraint.
T clear(T... args)
T end(T... args)
T iota(T... args)
T max(T... args)
contains data mapping from points to meshes.
constexpr auto get(span< E, S > s) -> decltype(s[N])
Definition span.hpp:602
T push_back(T... args)
T reserve(T... args)
T stable_sort(T... args)
int dataDims
The dimensionality of the data.
Definition Sample.hpp:45
Eigen::VectorXd values
Definition Sample.hpp:49