preCICE v3.1.2
Loading...
Searching...
No Matches
Mapping.cpp
Go to the documentation of this file.
1#include "Mapping.hpp"
2#include <boost/config.hpp>
3#include <ostream>
5#include "mesh/Utils.hpp"
6#include "utils/IntraComm.hpp"
7#include "utils/assertion.hpp"
8
9namespace precice::mapping {
10
12 Constraint constraint,
13 int dimensions,
14 bool requiresGradientData,
15 InitialGuessRequirement mappingType)
16 : _requiresGradientData(requiresGradientData),
17 _constraint(constraint),
18 _inputRequirement(MeshRequirement::UNDEFINED),
19 _outputRequirement(MeshRequirement::UNDEFINED),
20 _input(),
21 _output(),
22 _dimensions(dimensions),
23 _initialGuessRequirement(mappingType)
24{
25}
26
28 const mesh::PtrMesh &input,
29 const mesh::PtrMesh &output)
30{
31 _input = input;
33}
34
36{
37 return _input;
38}
39
41{
42 return _output;
43}
44
49
54
56{
57 PRECICE_ASSERT(requiresInitialGuess(), "This mapping isn't iterative, so it cannot have an initial guess.");
58 PRECICE_ASSERT(_initialGuess != nullptr, "The last solution wasn't provided.");
59 return _initialGuess->size() > 0;
60}
61
62const Eigen::VectorXd &Mapping::initialGuess() const
63{
64 PRECICE_ASSERT(requiresInitialGuess(), "This mapping isn't iterative, so it doesn't have an initial guess.");
65 PRECICE_ASSERT(_initialGuess != nullptr, "The last solution wasn't provided.");
66 return *_initialGuess;
67}
68
69Eigen::VectorXd &Mapping::initialGuess()
70{
71 PRECICE_ASSERT(requiresInitialGuess(), "This mapping isn't iterative, so it doesn't have an initial guess.");
72 PRECICE_ASSERT(_initialGuess != nullptr, "The last solution wasn't provided.");
73 return *_initialGuess;
74}
75
80
85
87{
88 return _input;
89}
90
92{
93 return _output;
94}
95
97 MeshRequirement requirement)
98{
99 _inputRequirement = requirement;
100}
101
103 MeshRequirement requirement)
104{
105 _outputRequirement = requirement;
106}
107
109{
110 return _dimensions;
111}
112
114{
116}
117
118void Mapping::map(int inputDataID, int outputDataID, Eigen::VectorXd &initialGuess)
119{
120 PRECICE_ASSERT(_initialGuess == nullptr);
122 map(inputDataID, outputDataID);
123 _initialGuess = nullptr;
124}
125
126void Mapping::map(int inputDataID,
127 int outputDataID)
128{
130 PRECICE_ASSERT(!requiresInitialGuess() || _initialGuess != nullptr, "Call the map version with lastSolution");
131
136 PRECICE_ASSERT(input()->data(inputDataID)->getDimensions() == output()->data(outputDataID)->getDimensions(),
137 input()->data(inputDataID)->getDimensions(), output()->data(outputDataID)->getDimensions());
138 PRECICE_ASSERT(input()->data(inputDataID)->values().size() / input()->data(inputDataID)->getDimensions() == static_cast<int>(input()->nVertices()),
139 input()->data(inputDataID)->values().size(), input()->data(inputDataID)->getDimensions(), input()->nVertices());
140 PRECICE_ASSERT(output()->data(outputDataID)->values().size() / output()->data(outputDataID)->getDimensions() == static_cast<int>(output()->nVertices()),
141 output()->data(outputDataID)->values().size(), output()->data(outputDataID)->getDimensions(), output()->nVertices());
142
143 time::Sample sample{input()->data(inputDataID)->getDimensions(),
144 input()->data(inputDataID)->values(),
145 input()->data(inputDataID)->gradients()};
146 map(sample, output()->data(outputDataID)->values());
147}
148
149void Mapping::map(const time::Sample &input, Eigen::VectorXd &output, Eigen::VectorXd &lastSolution)
150{
151 PRECICE_ASSERT(_initialGuess == nullptr);
152 _initialGuess = &lastSolution;
153 map(input, output);
154 _initialGuess = nullptr;
155}
156
157void Mapping::map(const time::Sample &input, Eigen::VectorXd &output)
158{
160 PRECICE_ASSERT(!requiresInitialGuess() || _initialGuess != nullptr, "Call the map version with lastSolution");
161
164 } else if (hasConstraint(CONSISTENT)) {
166 } else if (isScaledConsistent()) {
169 } else {
170 PRECICE_UNREACHABLE("Unknown mapping constraint.")
171 }
172}
173
174void Mapping::scaleConsistentMapping(const Eigen::VectorXd &input, Eigen::VectorXd &output, Mapping::Constraint constraint) const
175{
177
178 if (input.size() == 0 || output.size() == 0) {
179 return;
180 }
181 PRECICE_ASSERT(input.size() > 0 && output.size() > 0);
182
183 bool volumeMode = hasConstraint(SCALED_CONSISTENT_VOLUME);
184 logging::Logger _log{"mapping::Mapping"};
185 // Only serial participant is supported for scale-consistent mapping
187
188 // If rank is not empty and do not contain connectivity information, raise error
189 int spaceDimension = this->input()->getDimensions();
190 bool requiresEdges = (spaceDimension == 2 and !volumeMode);
191 bool requiresTriangles = (spaceDimension == 2 and volumeMode) or (spaceDimension == 3 and !volumeMode);
192 bool requiresTetra = (spaceDimension == 3 and volumeMode);
193
194 for (mesh::PtrMesh mesh : {this->input(), this->output()}) {
195 if (not mesh->empty()) {
196
197 PRECICE_CHECK(!(requiresEdges && mesh->edges().empty()), "Edges connectivity information is missing for the mesh \"{}\". "
198 "Scaled consistent mapping requires connectivity information.",
199 mesh->getName());
200
201 PRECICE_CHECK(!(requiresTriangles && mesh->triangles().empty()), "Triangles connectivity information is missing for the mesh \"{}\". "
202 "Scaled consistent mapping requires connectivity information.",
203 mesh->getName());
204
205 PRECICE_CHECK(!(requiresTetra && mesh->tetrahedra().empty()), "Tetrahedra connectivity information is missing for the mesh \"{}\". "
206 "Scaled consistent mapping requires connectivity information.",
207 mesh->getName());
208 }
209 }
210
211 const int valueDimensions = input.size() / this->input()->nVertices();
212
213 Eigen::VectorXd integralInput;
214 Eigen::VectorXd integralOutput;
215
216 // Integral is calculated on each direction separately
217 if (!volumeMode) {
218 integralInput = mesh::integrateSurface(this->input(), input);
219 integralOutput = mesh::integrateSurface(this->output(), output);
220 } else {
221 integralInput = mesh::integrateVolume(this->input(), input);
222 integralOutput = mesh::integrateVolume(this->output(), output);
223 }
224
225 // Create reshape the output values vector to matrix
226 Eigen::Map<Eigen::MatrixXd> outputValuesMatrix(output.data(), valueDimensions, output.size() / valueDimensions);
227
228 // Scale in each direction
229 // We cannot handle the case with zero output data and non-zero input data.
230 // To fulfill the constraint, we would need to scale the output data in such a way, that the integral sum of the input is preserved.
231 // That's not possible using a constant scaling factor as the output sum will always be zero. Here, we return 1 and emit a warning afterwards.
232 const Eigen::VectorXd scalingFactor = integralInput.binaryExpr(integralOutput, [](double lhs, double rhs) { return (rhs == 0.0) ? 1.0 : (lhs / rhs); });
233 PRECICE_DEBUG("Scaling factor in scaled-consistent mapping: {}", scalingFactor);
234
235 PRECICE_DEBUG("Scaling factor in scaled-consistent mapping: {}", scalingFactor);
236 outputValuesMatrix.array().colwise() *= scalingFactor.array();
237
238 // check whether the constraint is fulfilled
239 for (Eigen::Index i = 0; i < scalingFactor.size(); ++i) {
240 double consistency = scalingFactor[i] * integralOutput[i] - integralInput[i];
242 math::greater(std::abs(consistency), 0.0),
243 "Failed to fulfill consistency constraint of component {} for scaled-consistent mapping from mesh \"{}\" to mesh \"{}\". Consistency difference between input and scaled output is \"{}\".", i, this->input()->getName(), this->output()->getName(), consistency);
244 }
245} // namespace mapping
246
247bool Mapping::hasConstraint(const Constraint &constraint) const
248{
249 return (getConstraint() == constraint);
250}
251
253{
254 return _hasComputedMapping;
255}
256
261
263{
264 switch (lhs) {
268 return rhs == Mapping::MeshRequirement::FULL;
270 return false;
271 };
272 BOOST_UNREACHABLE_RETURN(false);
273}
274
276{
277 switch (val) {
279 out << "UNDEFINED";
280 break;
282 out << "VERTEX";
283 break;
285 out << "FULL";
286 break;
287 default:
288 PRECICE_ASSERT(false, "Implementation does not cover all cases");
289 };
290 return out;
291}
292
293} // namespace precice::mapping
std::ostream & out
#define PRECICE_WARN_IF(condition,...)
Definition LogMacros.hpp:21
#define PRECICE_DEBUG(...)
Definition LogMacros.hpp:64
#define PRECICE_CHECK(check,...)
Definition LogMacros.hpp:35
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
#define PRECICE_UNREACHABLE(...)
Definition assertion.hpp:95
This class provides a lightweight logger.
Definition Logger.hpp:16
bool requiresInitialGuess() const
Return true if the mapping requires an initial guess.
Definition Mapping.cpp:50
mesh::PtrMesh output() const
Returns pointer to output mesh.
Definition Mapping.cpp:91
Eigen::VectorXd & initialGuess()
Definition Mapping.cpp:69
Constraint _constraint
Determines whether mapping is consistent or conservative.
Definition Mapping.hpp:242
mesh::PtrMesh _input
Pointer to input mesh.
Definition Mapping.hpp:251
virtual void mapConsistent(const time::Sample &input, Eigen::VectorXd &output)=0
Maps data using a consistent constraint.
MeshRequirement
Specifies requirements for the input and output meshes of a mapping.
Definition Mapping.hpp:45
Constraint
Specifies additional constraints for a mapping.
Definition Mapping.hpp:29
Eigen::VectorXd * _initialGuess
Pointer to the initialGuess set and unset by map.
Definition Mapping.hpp:262
InitialGuessRequirement _initialGuessRequirement
The InitialGuessRequirement of the Mapping.
Definition Mapping.hpp:259
void setMeshes(const mesh::PtrMesh &input, const mesh::PtrMesh &output)
Sets input and output meshes carrying data to be mapped.
Definition Mapping.cpp:27
Mapping(Constraint constraint, int dimensions, bool requiresGradientData, InitialGuessRequirement initialGuessRequirement)
Constructor, takes mapping constraint.
Definition Mapping.cpp:11
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
bool isScaledConsistent() const
Returns true if mapping is a form of scaled consistent mapping.
Definition Mapping.cpp:257
virtual void scaleConsistentMapping(const Eigen::VectorXd &input, Eigen::VectorXd &output, Constraint type) const
Scales the consistently mapped output data such that the surface integral of the values on input mesh...
Definition Mapping.cpp:174
void setInputRequirement(MeshRequirement requirement)
Sets the mesh requirement for the input mesh.
Definition Mapping.cpp:96
bool hasInitialGuess() const
True if initialGuess().size() == 0.
Definition Mapping.cpp:55
mesh::PtrMesh _output
Pointer to output mesh.
Definition Mapping.hpp:254
bool requiresGradientData() const
Returns whether the mapping requires gradient data.
Definition Mapping.cpp:113
MeshRequirement _inputRequirement
Requirement on input mesh.
Definition Mapping.hpp:245
virtual std::string getName() const =0
Returns the name of the mapping method for logging purpose.
virtual void mapConservative(const time::Sample &input, Eigen::VectorXd &output)=0
Maps data using a conservative constraint.
MeshRequirement _outputRequirement
Requirement on output mesh.
Definition Mapping.hpp:248
const mesh::PtrMesh & getOutputMesh() const
Definition Mapping.cpp:40
void setOutputRequirement(MeshRequirement requirement)
Sets the mesh requirement for the output mesh.
Definition Mapping.cpp:102
const Eigen::VectorXd & initialGuess() const
Return the provided initial guess of a mapping using an initialGuess.
Definition Mapping.cpp:62
bool hasComputedMapping() const
Returns true, if the mapping has been computed.
Definition Mapping.cpp:252
void map(int inputDataID, int outputDataID)
Definition Mapping.cpp:126
MeshRequirement getInputRequirement() const
Returns the requirement on the input mesh.
Definition Mapping.cpp:76
const mesh::PtrMesh & getInputMesh() const
Definition Mapping.cpp:35
MeshRequirement getOutputRequirement() const
Returns the requirement on the output mesh.
Definition Mapping.cpp:81
virtual bool hasConstraint(const Constraint &constraint) const
Checks whether the mapping has the given constraint or not.
Definition Mapping.cpp:247
InitialGuessRequirement
Specifies whether the mapping requires an initial guess.
Definition Mapping.hpp:63
bool _requiresGradientData
Flag if gradient data is required for the mapping.
Definition Mapping.hpp:211
static bool isPrimary()
True if this process is running the primary rank.
Definition IntraComm.cpp:52
static bool isSecondary()
True if this process is running a secondary rank.
Definition IntraComm.cpp:57
contains data mapping from points to meshes.
bool operator<(Mapping::MeshRequirement lhs, Mapping::MeshRequirement rhs)
Definition Mapping.cpp:262
std::ostream & operator<<(std::ostream &out, Mapping::MeshRequirement val)
Definition Mapping.cpp:275
std::enable_if< std::is_arithmetic< Scalar >::value, bool >::type greater(Scalar A, Scalar B, Scalar tolerance=NUMERICAL_ZERO_DIFFERENCE)
Eigen::VectorXd integrateSurface(const PtrMesh &mesh, const Eigen::VectorXd &input)
Given the data and the mesh, this function returns the surface integral. Assumes no overlap exists fo...
Definition Utils.cpp:11
Eigen::VectorXd integrateVolume(const PtrMesh &mesh, const Eigen::VectorXd &input)
Given the data and the mesh, this function returns the volume integral. Assumes no overlap exists for...
Definition Utils.cpp:41
static precice::logging::Logger _log("precicec")