preCICE v3.1.1
Loading...
Searching...
No Matches
Utils.cpp
Go to the documentation of this file.
1#include <Eigen/Core>
2#include <mesh/Edge.hpp>
3#include <mesh/Mesh.hpp>
4#include <mesh/Utils.hpp>
5#include <utils/IntraComm.hpp>
6#include "utils/assertion.hpp"
7
8namespace precice::mesh {
9
11Eigen::VectorXd integrateSurface(const PtrMesh &mesh, const Eigen::VectorXd &input)
12{
13 PRECICE_ASSERT(mesh->nVertices() > 0);
14 const int meshDimensions = mesh->getDimensions();
15 const int valueDimensions = input.size() / mesh->nVertices();
16 const auto & values = input;
17 Eigen::VectorXd integral = Eigen::VectorXd::Zero(valueDimensions);
18
19 if (meshDimensions == 2) {
20 for (const auto &edge : mesh->edges()) {
21 int vertex1 = edge.vertex(0).getID() * valueDimensions;
22 int vertex2 = edge.vertex(1).getID() * valueDimensions;
23 for (int dim = 0; dim < valueDimensions; ++dim) {
24 integral(dim) += 0.5 * edge.getLength() * (values(vertex1 + dim) + values(vertex2 + dim));
25 }
26 }
27 } else {
28 for (const auto &face : mesh->triangles()) {
29 int vertex1 = face.vertex(0).getID() * valueDimensions;
30 int vertex2 = face.vertex(1).getID() * valueDimensions;
31 int vertex3 = face.vertex(2).getID() * valueDimensions;
32
33 for (int dim = 0; dim < valueDimensions; ++dim) {
34 integral(dim) += (face.getArea() / 3.0) * (values(vertex1 + dim) + values(vertex2 + dim) + values(vertex3 + dim));
35 }
36 }
37 }
38 return integral;
39}
40
41Eigen::VectorXd integrateVolume(const PtrMesh &mesh, const Eigen::VectorXd &input)
42{
43 PRECICE_ASSERT(mesh->nVertices() > 0);
44 const int meshDimensions = mesh->getDimensions();
45 const int valueDimensions = input.size() / mesh->nVertices();
46 const auto & values = input;
47 Eigen::VectorXd integral = Eigen::VectorXd::Zero(valueDimensions);
48 if (meshDimensions == 2) {
49 for (const auto &face : mesh->triangles()) {
50 int vertex1 = face.vertex(0).getID() * valueDimensions;
51 int vertex2 = face.vertex(1).getID() * valueDimensions;
52 int vertex3 = face.vertex(2).getID() * valueDimensions;
53 for (int dim = 0; dim < valueDimensions; ++dim) {
54 integral(dim) += (face.getArea() / 3.0) * (values(vertex1 + dim) + values(vertex2 + dim) + values(vertex3 + dim));
55 }
56 }
57 } else {
58 for (const auto &tetra : mesh->tetrahedra()) {
59 int vertex1 = tetra.vertex(0).getID() * valueDimensions;
60 int vertex2 = tetra.vertex(1).getID() * valueDimensions;
61 int vertex3 = tetra.vertex(2).getID() * valueDimensions;
62 int vertex4 = tetra.vertex(3).getID() * valueDimensions;
63
64 for (int dim = 0; dim < valueDimensions; ++dim) {
65 integral(dim) += (tetra.getVolume() / 4.0) * (values(vertex1 + dim) + values(vertex2 + dim) + values(vertex3 + dim) + values(vertex4 + dim));
66 }
67 }
68 }
69 return integral;
70}
71
72} // namespace precice::mesh
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
int getDimensions() const
Definition Mesh.cpp:98
std::size_t nVertices() const
Returns the number of vertices.
Definition Mesh.cpp:63
TetraContainer & tetrahedra()
Returns modifiable container holding all tetrahedra.
Definition Mesh.cpp:93
TriangleContainer & triangles()
Returns modifiable container holding all triangles.
Definition Mesh.cpp:78
EdgeContainer & edges()
Returns modifiable container holding all edges.
Definition Mesh.cpp:68
provides Mesh, Data and primitives.
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