preCICE v3.1.2
Loading...
Searching...
No Matches
LinearCellInterpolationMappingTest.cpp
Go to the documentation of this file.
1#include <Eigen/Core>
2#include <algorithm>
3#include <memory>
4#include "logging/LogMacros.hpp"
5#include "mapping/LinearCellInterpolationMapping.hpp"
6#include "mapping/Mapping.hpp"
7#include "math/constants.hpp"
8#include "mesh/Data.hpp"
9#include "mesh/Mesh.hpp"
10#include "mesh/SharedPointer.hpp"
11#include "mesh/Utils.hpp"
12#include "mesh/Vertex.hpp"
13#include "testing/TestContext.hpp"
14#include "testing/Testing.hpp"
15
16using namespace precice;
17using namespace precice::mesh;
18
19BOOST_AUTO_TEST_SUITE(MappingTests)
20BOOST_AUTO_TEST_SUITE(LinearCellInterpolationMapping)
21
22BOOST_AUTO_TEST_CASE(Consistent)
23{
24 PRECICE_TEST(1_rank);
25 int dimensions = 2;
26 using testing::equals;
27
28 PtrMesh inMesh(new Mesh("InMesh", dimensions, testing::nextMeshID()));
29 PtrData inDataScalar = inMesh->createData("InDataScalar", 1, 0_dataID);
30 int inDataScalarID = inDataScalar->getID();
31
32 Vertex &inVertexA = inMesh->createVertex(Eigen::Vector2d(0.0, 0.0));
33 Vertex &inVertexB = inMesh->createVertex(Eigen::Vector2d(1.0, 0.0));
34 Vertex &inVertexC = inMesh->createVertex(Eigen::Vector2d(0.0, 1.0));
35
36 inMesh->allocateDataValues();
37
38 Edge &inEdge0 = inMesh->createEdge(inVertexA, inVertexB);
39 Edge &inEdge1 = inMesh->createEdge(inVertexB, inVertexC);
40 Edge &inEdge2 = inMesh->createEdge(inVertexC, inVertexA);
41
42 inMesh->createTriangle(inEdge0, inEdge1, inEdge2);
43 Eigen::VectorXd &inValuesScalar = inDataScalar->values();
44 inValuesScalar << 1.0, 2.0, 3.0;
45
46 BOOST_CHECK(!inMesh->edges().empty());
47 /*
48 Should pass: triangles are discarded when the mapping doesn't require the in Participant::setMeshTriangle,
49 but Mesh::createTriangle doesn't check it. This is not an integration test!
50 */
51
52 BOOST_CHECK(!inMesh->triangles().empty());
53 PtrMesh outMesh(new Mesh("OutMesh", dimensions, testing::nextMeshID()));
54 PtrData outDataScalar = outMesh->createData("OutDataScalar", 1, 2_dataID);
55 int outDataScalarID = outDataScalar->getID();
56
57 // All vertices to test
58 // Center of triangle = average
59 outMesh->createVertex(Eigen::Vector2d::Constant(1.0 / 3.0));
60 // Exact mapping if grid is matching
61 outMesh->createVertex(Eigen::Vector2d(0.0, 0.0));
62 outMesh->createVertex(Eigen::Vector2d(1.0, 0.0));
63 outMesh->createVertex(Eigen::Vector2d(0.0, 1.0));
64 // Fallback on NP when slightly outside
65 // AB: exact middle (slightly outside and on side of A), BC: 2/3 on side B. CA: check fall-back on edge if out of domain
66 outMesh->createVertex(Eigen::Vector2d(0.49, -0.01));
67 outMesh->createVertex(Eigen::Vector2d(2.5 / 3, 1. / 3));
68 outMesh->createVertex(Eigen::Vector2d(-10.0, 0.25));
69 // Check fall back on nearest neighbor
70 outMesh->createVertex(Eigen::Vector2d(-0.1, -0.1)); // Fall back to NN expected
71 outMesh->createVertex(Eigen::Vector2d(2.5, -1.0));
72 outMesh->createVertex(Eigen::Vector2d(2.5, 10.0));
73
74 outMesh->allocateDataValues();
75
76 // Setup mapping with mapping coordinates and geometry used
77 precice::mapping::LinearCellInterpolationMapping mapping(mapping::Mapping::CONSISTENT, dimensions);
78 mapping.setMeshes(inMesh, outMesh);
79 BOOST_TEST(mapping.hasComputedMapping() == false);
80 mapping.computeMapping();
81 mapping.map(inDataScalarID, outDataScalarID);
82 const Eigen::VectorXd &outValuesScalar = outDataScalar->values();
83 BOOST_TEST(mapping.hasComputedMapping() == true);
84
85 // Check expected
86 Eigen::VectorXd expected(outMesh->nVertices());
87 expected << 2.0, 1.0, 2.0, 3.0, 1.49, 2.25, 1.5, 1.0, 2.0, 3.0;
88 BOOST_CHECK(equals(expected, outValuesScalar));
89}
90
91BOOST_AUTO_TEST_CASE(Conservative)
92{
93 PRECICE_TEST(1_rank);
94 int dimensions = 2;
95 using testing::equals;
96
97 /* Send forces on some points, map to a triangle and check correct distribution.
98 We only apply forces inside the triangle or at its boundary.
99 Behavior when projecting is not ideal, see #1304
100 */
101
102 // Map to a triangle ABC which represents the lower-left half of a unit square
103 const double forceOnMid = 1.0;
104 const double forceOnMidAB = 2.0;
105 const double forceOnA = 10.0;
106 const double forceOnC = 7.0;
107 const double unbalancedforceOnBC = 3.0; // 75% on B, 25% on C
108 const double netForce = forceOnMid + forceOnMidAB + forceOnA + forceOnC + unbalancedforceOnBC;
109
110 PtrMesh inMesh(new Mesh("InMesh", dimensions, testing::nextMeshID()));
111 PtrData inDataScalar = inMesh->createData("InDataScalar", 1, 0_dataID);
112 int inDataScalarID = inDataScalar->getID();
113
114 inMesh->createVertex(Eigen::Vector2d(1. / 3, 1. / 3)); // Mid
115 inMesh->createVertex(Eigen::Vector2d(0.5, 0.0)); // Mid AB
116 inMesh->createVertex(Eigen::Vector2d(0.0, 0.0)); // A
117 inMesh->createVertex(Eigen::Vector2d(0.0, 1.0)); // C
118 inMesh->createVertex(Eigen::Vector2d(0.75, 0.25)); // Along BC
119
120 inMesh->allocateDataValues();
121 Eigen::VectorXd &inValuesScalar = inDataScalar->values();
122 inValuesScalar << forceOnMid, forceOnMidAB, forceOnA, forceOnC, unbalancedforceOnBC;
123
124 // Create output mesh (the triangle ABC)
125 PtrMesh outMesh(new Mesh("OutMesh", dimensions, testing::nextMeshID()));
126 PtrData outDataScalar = outMesh->createData("OutDataScalar", 1, 2_dataID);
127 int outDataScalarID = outDataScalar->getID();
128
129 Vertex &outVertexA = outMesh->createVertex(Eigen::Vector2d(0.0, 0.0));
130 Vertex &outVertexB = outMesh->createVertex(Eigen::Vector2d(1.0, 0.0));
131 Vertex &outVertexC = outMesh->createVertex(Eigen::Vector2d(0.0, 1.0));
132
133 Edge &outEdge0 = outMesh->createEdge(outVertexA, outVertexB);
134 Edge &outEdge1 = outMesh->createEdge(outVertexB, outVertexC);
135 Edge &outEdge2 = outMesh->createEdge(outVertexC, outVertexA);
136
137 outMesh->createTriangle(outEdge0, outEdge1, outEdge2);
138 BOOST_CHECK(!outMesh->edges().empty());
139 BOOST_CHECK(!outMesh->triangles().empty());
140
141 outMesh->allocateDataValues();
142
143 // Setup mapping with mapping coordinates and geometry used
144 precice::mapping::LinearCellInterpolationMapping mapping(mapping::Mapping::CONSERVATIVE, dimensions);
145 mapping.setMeshes(inMesh, outMesh);
146 BOOST_TEST(mapping.hasComputedMapping() == false);
147 mapping.computeMapping();
148 mapping.map(inDataScalarID, outDataScalarID);
149 const Eigen::VectorXd &outValuesScalar = outDataScalar->values();
150 BOOST_TEST(mapping.hasComputedMapping() == true);
151
152 // Check expected
153 Eigen::VectorXd expected(outMesh->nVertices());
154 const double expectedA = forceOnMid / 3 + forceOnMidAB / 2 + forceOnA;
155 const double expectedB = forceOnMid / 3 + forceOnMidAB / 2 + unbalancedforceOnBC * 0.75;
156 const double expectedC = forceOnMid / 3 + forceOnC + unbalancedforceOnBC * 0.25;
157 expected << expectedA, expectedB, expectedC;
158 // Check expected value, and conservation
159 BOOST_CHECK(equals(expected, outValuesScalar));
160 BOOST_CHECK(equals(netForce, outValuesScalar.sum()));
161}
162
163BOOST_AUTO_TEST_CASE(ConsistentOneTetra3D)
164{
165 PRECICE_TEST(1_rank);
166 int dimensions = 3;
167 using testing::equals;
168
169 PtrMesh inMesh(new Mesh("InMesh", dimensions, testing::nextMeshID()));
170 PtrData inDataScalar = inMesh->createData("InDataScalar", 1, 0_dataID);
171 int inDataScalarID = inDataScalar->getID();
172
173 Vertex &inVertexA = inMesh->createVertex(Eigen::Vector3d(0.0, 0.0, 0.0));
174 Vertex &inVertexB = inMesh->createVertex(Eigen::Vector3d(1.0, 0.0, 0.0));
175 Vertex &inVertexC = inMesh->createVertex(Eigen::Vector3d(0.0, 1.0, 0.0));
176 Vertex &inVertexD = inMesh->createVertex(Eigen::Vector3d(0.0, 0.0, 1.0));
177
178 inMesh->allocateDataValues();
179
180 // Create a tetra
181 inMesh->createTetrahedron(inVertexA, inVertexB, inVertexC, inVertexD);
182 // Create triangle in the plane z = 0
183 inMesh->createTriangle(inVertexA, inVertexB, inVertexC);
184 // Add edge BD to check fall-back on edge
185 inMesh->createEdge(inVertexB, inVertexD);
186
187 Eigen::VectorXd &inValuesScalar = inDataScalar->values();
188 inValuesScalar << 1.0, 2.0, 3.0, 4.0; //1 + x + 2y + 3z
189
190 BOOST_CHECK(!inMesh->tetrahedra().empty());
191 PtrMesh outMesh(new Mesh("OutMesh", dimensions, testing::nextMeshID()));
192 PtrData outDataScalar = outMesh->createData("OutDataScalar", 1, 2_dataID);
193 int outDataScalarID = outDataScalar->getID();
194
195 // Center and the 4 vertices (expected: average then 4 values)
196 outMesh->createVertex(Eigen::Vector3d::Constant(0.25));
197 outMesh->createVertex(Eigen::Vector3d(0.0, 0.0, 0.0));
198 outMesh->createVertex(Eigen::Vector3d(1.0, 0.0, 0.0));
199 outMesh->createVertex(Eigen::Vector3d(0.0, 1.0, 0.0));
200 outMesh->createVertex(Eigen::Vector3d(0.0, 0.0, 1.0));
201 // Point below the triangle ABC -> fallback to triangle. Expected projection on triangle.
202 outMesh->createVertex(Eigen::Vector3d(1.0 / 3, 1.0 / 3, -0.1));
203 // Point close to the triangle ACD (which isn't set!).
204 // Wanted behavior: NN => 3.0. Incorrect behavior: 3.6 in case of fall-back to edge. See issue #1304
205 outMesh->createVertex(Eigen::Vector3d(-0.1, 0.8, 0.5));
206 // Point inside the triangle BCD (not set) -> Check it is inside the tetra. Expected: 0.2*2+0.3*3+0.5*4 = 3.3
207 outMesh->createVertex(Eigen::Vector3d(0.2, 0.3, 0.5));
208 // Point on the the edge BD for fall-back. Expected: 0.4*3 + 0.6*4 = 3.6
209 outMesh->createVertex(Eigen::Vector3d(0, 0.4, 0.6));
210
211 outMesh->allocateDataValues();
212
213 // Setup mapping with mapping coordinates and geometry used
214 precice::mapping::LinearCellInterpolationMapping mapping(mapping::Mapping::CONSISTENT, dimensions);
215 mapping.setMeshes(inMesh, outMesh);
216 BOOST_TEST(mapping.hasComputedMapping() == false);
217 mapping.computeMapping();
218 mapping.map(inDataScalarID, outDataScalarID);
219 const Eigen::VectorXd &outValuesScalar = outDataScalar->values();
220 BOOST_TEST(mapping.hasComputedMapping() == true);
221
222 // Check expected
223 Eigen::VectorXd expected(outMesh->nVertices());
224 expected << 2.5, 1.0, 2.0, 3.0, 4.0, 2.0, 3.0, 3.3, 3.6;
225 BOOST_CHECK(equals(expected, outValuesScalar));
226}
227
228BOOST_AUTO_TEST_CASE(ConservativeOneTetra3D)
229{
230 PRECICE_TEST(1_rank);
231 int dimensions = 3;
232 using testing::equals;
233
234 /* Send forces on some points, map to a tetra and check correct distribution.
235 We only apply forces inside the tetra or at its boundary.
236 Behavior when projecting is not ideal, see #1304
237 */
238
239 const double forceOnMid = 1.0;
240 const double forceOnMidAB = 2.0;
241 const double forceOnA = 10.0;
242 const double forceOnC = 7.0;
243 const double unbalancedforceOnBC = 3.0; // 75% on B, 25% on C
244 const double unbalancedInternalForce = 11.0; // With barycentric coordinates (0.4, 0.1, 0.2, 0.3)
245 const double netForce = forceOnMid + forceOnMidAB + forceOnA + forceOnC + unbalancedforceOnBC + unbalancedInternalForce;
246
247 PtrMesh inMesh(new Mesh("InMesh", dimensions, testing::nextMeshID()));
248 PtrData inDataScalar = inMesh->createData("InDataScalar", 1, 0_dataID);
249 int inDataScalarID = inDataScalar->getID();
250
251 inMesh->createVertex(Eigen::Vector3d(0.25, 0.25, 0.25)); // Mid
252 inMesh->createVertex(Eigen::Vector3d(0.5, 0.0, 0.0)); // Mid AB
253 inMesh->createVertex(Eigen::Vector3d(0.0, 0.0, 0.0)); // A
254 inMesh->createVertex(Eigen::Vector3d(0.0, 1.0, 0.0)); // C
255 inMesh->createVertex(Eigen::Vector3d(0.75, 0.25, 0.0)); // Along BC
256 inMesh->createVertex(Eigen::Vector3d(0.1, 0.2, 0.3));
257
258 inMesh->allocateDataValues();
259 Eigen::VectorXd &inValuesScalar = inDataScalar->values();
260 inValuesScalar << forceOnMid, forceOnMidAB, forceOnA, forceOnC, unbalancedforceOnBC, unbalancedInternalForce;
261
262 // Create output mesh (the tetra ABCD)
263 PtrMesh outMesh(new Mesh("OutMesh", dimensions, testing::nextMeshID()));
264 PtrData outDataScalar = outMesh->createData("OutDataScalar", 1, 2_dataID);
265 int outDataScalarID = outDataScalar->getID();
266
267 Vertex &outVertexA = outMesh->createVertex(Eigen::Vector3d(0.0, 0.0, 0.0));
268 Vertex &outVertexB = outMesh->createVertex(Eigen::Vector3d(1.0, 0.0, 0.0));
269 Vertex &outVertexC = outMesh->createVertex(Eigen::Vector3d(0.0, 1.0, 0.0));
270 Vertex &outVertexD = outMesh->createVertex(Eigen::Vector3d(0.0, 0.0, 1.0));
271
272 outMesh->createTetrahedron(outVertexA, outVertexB, outVertexC, outVertexD);
273 BOOST_CHECK(!outMesh->tetrahedra().empty());
274
275 outMesh->allocateDataValues();
276
277 // Setup mapping with mapping coordinates and geometry used
278 precice::mapping::LinearCellInterpolationMapping mapping(mapping::Mapping::CONSERVATIVE, dimensions);
279 mapping.setMeshes(inMesh, outMesh);
280 BOOST_TEST(mapping.hasComputedMapping() == false);
281 mapping.computeMapping();
282 mapping.map(inDataScalarID, outDataScalarID);
283 const Eigen::VectorXd &outValuesScalar = outDataScalar->values();
284 BOOST_TEST(mapping.hasComputedMapping() == true);
285
286 // Check expected
287 Eigen::VectorXd expected(outMesh->nVertices());
288 const double expectedA = forceOnMid / 4 + forceOnMidAB / 2 + forceOnA + 0.4 * unbalancedInternalForce;
289 const double expectedB = forceOnMid / 4 + forceOnMidAB / 2 + unbalancedforceOnBC * 0.75 + 0.1 * unbalancedInternalForce;
290 const double expectedC = forceOnMid / 4 + forceOnC + unbalancedforceOnBC * 0.25 + 0.2 * unbalancedInternalForce;
291 const double expectedD = forceOnMid / 4 + 0.3 * unbalancedInternalForce;
292
293 expected << expectedA, expectedB, expectedC, expectedD;
294 // Check expected value, and conservation
295 BOOST_CHECK(equals(expected, outValuesScalar));
296 BOOST_CHECK(equals(netForce, outValuesScalar.sum()));
297}
298
299BOOST_AUTO_TEST_SUITE_END()
300BOOST_AUTO_TEST_SUITE_END()
BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(Consistent)
Definition LinearCellInterpolationMappingTest.cpp:22
BOOST_AUTO_TEST_SUITE
BOOST_AUTO_TEST_SUITE(PreProcess)
BOOST_AUTO_TEST_SUITE_END
BOOST_AUTO_TEST_SUITE_END()
PRECICE_TEST
#define PRECICE_TEST(...)
Definition Testing.hpp:27
precice::mapping::LinearCellInterpolationMapping
Mapping using orthogonal projection to nearest triangle/edge/vertex and linear interpolation from pro...
Definition LinearCellInterpolationMapping.hpp:13
precice::mapping::LinearCellInterpolationMapping::computeMapping
void computeMapping() final override
Computes the projections and interpolation relations.
Definition LinearCellInterpolationMapping.cpp:31
precice::mapping::Mapping::setMeshes
void setMeshes(const mesh::PtrMesh &input, const mesh::PtrMesh &output)
Sets input and output meshes carrying data to be mapped.
Definition Mapping.cpp:27
precice::mapping::Mapping::hasComputedMapping
bool hasComputedMapping() const
Returns true, if the mapping has been computed.
Definition Mapping.cpp:252
precice::mapping::Mapping::map
void map(int inputDataID, int outputDataID)
Definition Mapping.cpp:126
precice::mesh::Data::getID
DataID getID() const
Returns the ID of the data set (supposed to be unique).
Definition Data.cpp:93
precice::mesh::Data::values
Eigen::VectorXd & values()
Returns a reference to the data values.
Definition Data.cpp:28
precice::mesh::Edge
Linear edge of a mesh, defined by two Vertex objects.
Definition Edge.hpp:16
precice::mesh::Mesh
Container and creator for meshes.
Definition Mesh.hpp:39
precice::mesh::Mesh::createTriangle
Triangle & createTriangle(Edge &edgeOne, Edge &edgeTwo, Edge &edgeThree)
Creates and initializes a Triangle object.
Definition Mesh.cpp:119
precice::mesh::Mesh::nVertices
std::size_t nVertices() const
Returns the number of vertices.
Definition Mesh.cpp:63
precice::mesh::Mesh::tetrahedra
TetraContainer & tetrahedra()
Returns modifiable container holding all tetrahedra.
Definition Mesh.cpp:93
precice::mesh::Mesh::createData
PtrData & createData(const std::string &name, int dimension, DataID id, int waveformDegree=time::Time::DEFAULT_WAVEFORM_DEGREE)
Create only data for vertex.
Definition Mesh.cpp:151
precice::mesh::Mesh::createTetrahedron
Tetrahedron & createTetrahedron(Vertex &vertexOne, Vertex &vertexTwo, Vertex &vertexThree, Vertex &vertexFour)
Creates and initializes a Tetrahedron object.
Definition Mesh.cpp:141
precice::mesh::Mesh::triangles
TriangleContainer & triangles()
Returns modifiable container holding all triangles.
Definition Mesh.cpp:78
precice::mesh::Mesh::createEdge
Edge & createEdge(Vertex &vertexOne, Vertex &vertexTwo)
Creates and initializes an Edge object.
Definition Mesh.cpp:111
precice::mesh::Mesh::edges
EdgeContainer & edges()
Returns modifiable container holding all edges.
Definition Mesh.cpp:68
precice::mesh::Mesh::createVertex
Vertex & createVertex(const Eigen::VectorXd &coords)
Creates and initializes a Vertex object.
Definition Mesh.cpp:103
precice::mesh::Mesh::allocateDataValues
void allocateDataValues()
Allocates memory for the vertex data values and corresponding gradient values.
Definition Mesh.cpp:233
precice::mesh::Vertex
Vertex of a mesh.
Definition Vertex.hpp:16
std::deque::empty
T empty(T... args)
precice::mesh
provides Mesh, Data and primitives.
Definition ScaleActionTest.cpp:20
precice::testing::equals
boost::test_tools::predicate_result equals(const std::vector< float > &VectorA, const std::vector< float > &VectorB, float tolerance)
equals to be used in tests. Compares two std::vectors using a given tolerance. Prints both operands o...
Definition Testing.cpp:65
precice
Main namespace of the precice library.
Definition Acceleration.cpp:5