preCICE v3.1.1
Loading...
Searching...
No Matches
AxialGeoMultiscaleMapping.cpp
Go to the documentation of this file.
2#include <Eigen/Core>
3#include "mesh/Mesh.hpp"
4
5namespace precice::mapping {
6
8 Constraint constraint,
9 int dimensions,
10 MultiscaleType type,
11 MultiscaleAxis axis,
12 double radius)
13 : Mapping(constraint, dimensions, false, Mapping::InitialGuessRequirement::None),
14 _type(type),
15 _axis(axis),
16 _radius(radius)
17{
20}
21
23{
24 PRECICE_TRACE(output()->nVertices());
25
26 PRECICE_ASSERT(input().get() != nullptr);
27 PRECICE_ASSERT(output().get() != nullptr);
28
29 if (getConstraint() == CONSISTENT) {
30 PRECICE_DEBUG("Compute consistent mapping");
32 PRECICE_CHECK(input()->nVertices() == 1, "You can only define an axial geometric multiscale mapping of type spread from a mesh with exactly one vertex.");
33
34 /* When we add support for 1D meshes (https://github.com/precice/precice/issues/1669),
35 we should check for valid dimensions combination.
36
37 const int inDataDimensions = input()->getDimensions();
38 const int outDataDimensions = output()->getDimensions();
39
40 PRECICE_CHECK(input()->getDimensions() == 1, "The input mesh on an axial geometric multiscale mapping can only be 1D at the moment, but it was defined to be {}.", input()->getDimensions())
41 PRECICE_CHECK(output()->getDimensions() == 3, "The output mesh on an axial geometric multiscale mapping can only be 3D at the moment, but it was defined to be {}.", input()->getDimensions())
42 */
43 const int outDataDimensions = 3;
44
45 int effectiveCoordinate = static_cast<std::underlying_type_t<MultiscaleType>>(_axis); // Convert enum struct to int
46 PRECICE_ASSERT(effectiveCoordinate == static_cast<std::underlying_type_t<MultiscaleType>>(MultiscaleAxis::X) ||
47 effectiveCoordinate == static_cast<std::underlying_type_t<MultiscaleType>>(MultiscaleAxis::Y) ||
48 effectiveCoordinate == static_cast<std::underlying_type_t<MultiscaleType>>(MultiscaleAxis::Z),
49 "Unknown multiscale axis type.");
50
51 // compute distances between 1D vertex and 3D vertices
52 mesh::Vertex & v0 = input()->vertex(0);
53 size_t const outSize = output()->nVertices();
54 constexpr double distance_to_radius_threshold = 1.05;
55
57 _vertexDistances.reserve(output()->nVertices());
58
59 for (size_t i = 0; i < outSize; i++) {
60 Eigen::VectorXd difference(outDataDimensions);
61 difference = v0.getCoords();
62 difference -= output()->vertex(i).getCoords();
63 double distance_to_radius = difference.norm() / _radius;
64 PRECICE_CHECK(distance_to_radius <= distance_to_radius_threshold, "Output mesh has vertices that do not coincide with the geometric multiscale interface defined by the input mesh. Ratio of vertex distance to radius is {} (which is larger than the assumed threshold of distance_to_radius_threshold).", distance_to_radius);
65 _vertexDistances.push_back(distance_to_radius);
66 }
67 } else {
69 PRECICE_CHECK(output()->nVertices() == 1, "You can only define an axial geometric multiscale mapping of type spread from a mesh with exactly one vertex.");
70 // Nothing to do here: A consistent collect mapping only averages all the values, independently of their locations, and this is done in the mapConsistent() method.
71 }
72 } else {
74 PRECICE_UNREACHABLE("Axial conservative geometric multiscale mapping is not implemented");
75 PRECICE_DEBUG("Compute conservative mapping");
76 }
78}
79
86
87void AxialGeoMultiscaleMapping::mapConservative(const time::Sample &inData, Eigen::VectorXd &outData)
88{
90 PRECICE_UNREACHABLE("Axial conservative geometric multiscale mapping is not implemented");
91 PRECICE_DEBUG("Map conservative");
92}
93
94void AxialGeoMultiscaleMapping::mapConsistent(const time::Sample &inData, Eigen::VectorXd &outData)
95{
97
98 const int inDataDimensions = inData.dataDims;
99 const Eigen::VectorXd &inputValues = inData.values;
100 Eigen::VectorXd & outputValues = outData;
101 // TODO: check if this needs to change when access to mesh dimension is possible
102 const int outDataDimensions = outData.size() / output()->nVertices();
103
104 int effectiveCoordinate = static_cast<std::underlying_type_t<MultiscaleType>>(_axis);
105 PRECICE_ASSERT(effectiveCoordinate == static_cast<std::underlying_type_t<MultiscaleType>>(MultiscaleAxis::X) ||
106 effectiveCoordinate == static_cast<std::underlying_type_t<MultiscaleType>>(MultiscaleAxis::Y) ||
107 effectiveCoordinate == static_cast<std::underlying_type_t<MultiscaleType>>(MultiscaleAxis::Z),
108 "Unknown multiscale axis type.");
109
110 // Check that the number of values for the input and output is right according to their dimensions
111 PRECICE_ASSERT((inputValues.size() / static_cast<std::size_t>(inDataDimensions) == input()->nVertices()),
112 inputValues.size(), inDataDimensions, input()->nVertices());
113 PRECICE_ASSERT((outputValues.size() / static_cast<std::size_t>(outDataDimensions) == output()->nVertices()),
114 outputValues.size(), outDataDimensions, output()->nVertices());
115
116 // We currently don't support 1D data, so we need that the user specifies data of the same dimensions on both sides
117 PRECICE_ASSERT(inDataDimensions == outDataDimensions);
118
119 PRECICE_DEBUG("Map consistent");
121 /*
122 3D vertices are assigned a value based on distance from the 1D vertex.
123 Currently, a Hagen-Poiseuille profile determines the velocity value.
124 */
125 PRECICE_ASSERT(input()->nVertices() == 1);
126 size_t const outSize = output()->nVertices();
127
128 for (size_t i = 0; i < outSize; i++) {
129 PRECICE_ASSERT(static_cast<size_t>((i * outDataDimensions) + effectiveCoordinate) < static_cast<size_t>(outputValues.size()), ((i * outDataDimensions) + effectiveCoordinate), outputValues.size());
130 // When adding support for 2D, remember that this should be 1.5 * inputValues(effectiveCoordinate) * (1 - (_vertexDistances[i] * _vertexDistances[i]));
131 outputValues((i * outDataDimensions) + effectiveCoordinate) = 2 * inputValues(effectiveCoordinate) * (1 - (_vertexDistances[i] * _vertexDistances[i]));
132 }
133 } else {
135 /*
136 1D vertex is assigned the averaged value over all 3D vertices at the outlet,
137 but only of the effectiveCoordinate component of the velocity vector.
138 */
139 PRECICE_ASSERT(output()->nVertices() == 1);
140 outputValues(effectiveCoordinate) = 0;
141 size_t const inSize = input()->nVertices();
142 for (size_t i = 0; i < inSize; i++) {
143 PRECICE_ASSERT(static_cast<size_t>((i * inDataDimensions) + effectiveCoordinate) < static_cast<size_t>(inputValues.size()), ((i * inDataDimensions) + effectiveCoordinate), inputValues.size());
144 outputValues(effectiveCoordinate) += inputValues((i * inDataDimensions) + effectiveCoordinate);
145 }
146 outputValues(effectiveCoordinate) = outputValues(effectiveCoordinate) / inSize;
147 }
148}
149
151{
153
155
156 if (getConstraint() == CONSISTENT) {
157 PRECICE_ASSERT(_type == MultiscaleType::SPREAD, "Not yet implemented");
158 PRECICE_ASSERT(input()->nVertices() == 1);
159
160 input()->tagAll();
161
162 } else {
164 PRECICE_UNREACHABLE("Axial conservative geometric multiscale mapping is not implemented");
165 }
166
167 clear();
168}
169
171{
173 // no operation needed here for the moment
174}
175
177{
178 return "axial-geomultiscale";
179}
180
181} // namespace precice::mapping
#define PRECICE_DEBUG(...)
Definition LogMacros.hpp:64
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
#define PRECICE_CHECK(check,...)
Definition LogMacros.hpp:35
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
#define PRECICE_UNREACHABLE(...)
Definition assertion.hpp:95
void mapConservative(const time::Sample &inData, Eigen::VectorXd &outData) override
Maps data using a conservative constraint.
double _radius
radius of the "tube" from or to which the data is mapped, i.e., radius of the circular interface betw...
MultiscaleType
Geometric multiscale nature of the mapping (spread or collect).
std::string getName() const final override
Returns name of the mapping.
AxialGeoMultiscaleMapping(Constraint constraint, int dimensions, MultiscaleType type, MultiscaleAxis axis, double radius)
Constructor.
MultiscaleAxis _axis
main axis along which axial geometric multiscale coupling happens
void tagMeshSecondRound() override
Method used by partition. Tags vertices that can be filtered out.
void tagMeshFirstRound() override
Method used by partition. Tags vertices that could be owned by this rank.
void computeMapping() override
Takes care of compute-heavy operations needed only once to set up the mapping.
void clear() override
Removes a computed mapping.
MultiscaleType _type
type of mapping, namely spread or collect
std::vector< double > _vertexDistances
computed vertex distances to map data from input vertex to output vertices
void mapConsistent(const time::Sample &inData, Eigen::VectorXd &outData) override
Maps data using a consistent constraint.
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
Vertex of a mesh.
Definition Vertex.hpp:16
Eigen::VectorXd getCoords() const
Returns the coordinates of the vertex.
Definition Vertex.hpp:116
T clear(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)
int dataDims
The dimensionality of the data.
Definition Sample.hpp:45
Eigen::VectorXd values
Definition Sample.hpp:49