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