preCICE v3.1.1
Loading...
Searching...
No Matches
WatchIntegral.cpp
Go to the documentation of this file.
1#include "WatchIntegral.hpp"
2#include <Eigen/Core>
4#include "mesh/Data.hpp"
5#include "mesh/Edge.hpp"
6#include "mesh/Mesh.hpp"
8#include "mesh/Triangle.hpp"
9#include "mesh/Utils.hpp"
10#include "mesh/Vertex.hpp"
11#include "utils/IntraComm.hpp"
12
13namespace precice::impl {
14
16 mesh::PtrMesh meshToWatch,
17 const std::string &exportFilename,
18 bool isScalingOn)
19 : _mesh(std::move(meshToWatch)),
20 _txtWriter(exportFilename),
21 _isScalingOn(isScalingOn)
22{
24
27 }
28
29 for (size_t i = 0; i < _mesh->data().size(); i++) {
30 _dataToExport.push_back(_mesh->data()[i]);
31
33 if (_dataToExport[i]->getDimensions() > 1) {
34
35 io::TXTTableWriter::DataType vectorType = _dataToExport[i]->getDimensions() == 2
38 _txtWriter.addData(_dataToExport[i]->getName(), vectorType);
39 } else {
41 }
42 }
43 }
44}
45
47{
48 // Do not add surface area column if there is no connectivity
49 if ((not utils::IntraComm::isSecondary()) and (not _mesh->edges().empty())) {
51 }
52 PRECICE_WARN_IF(_isScalingOn && _mesh->edges().empty(),
53 "Watch-integral is configured with scaling option on; "
54 "however, mesh {} does not contain connectivity information. "
55 "Therefore, the integral will be calculated without scaling.",
56 _mesh->getName());
57}
58
60 double time)
61{
62
64 _txtWriter.writeData("Time", time);
65 }
66
67 for (auto &elem : _dataToExport) {
68 const int dataDimensions = elem->getDimensions();
69 auto integral = calculateIntegral(elem);
70
71 if (utils::IntraComm::getSize() > 1) {
72 Eigen::VectorXd valueRecv = Eigen::VectorXd::Zero(dataDimensions);
73 utils::IntraComm::reduceSum(integral, valueRecv);
74 integral = std::move(valueRecv);
75 }
77 if (dataDimensions == 1) {
78 _txtWriter.writeData(elem->getName(), integral[0]);
79 } else if (dataDimensions == 2) {
80 _txtWriter.writeData(elem->getName(), Eigen::Vector2d(integral));
81 } else {
82 _txtWriter.writeData(elem->getName(), Eigen::Vector3d(integral));
83 }
84 }
85 }
86
87 // Calculate surface area only if there is connectivity information
88 if (not _mesh->edges().empty()) {
89 double surfaceArea = calculateSurfaceArea();
90 if (utils::IntraComm::getSize() > 1) {
91 double surfaceAreaSum = 0.0;
92 utils::IntraComm::reduceSum(surfaceArea, surfaceAreaSum);
93 surfaceArea = surfaceAreaSum;
94 }
96 _txtWriter.writeData("SurfaceArea", surfaceArea);
97 }
98 } else {
99 // Empty partitions should also call reduceSum to prevent deadlock
100 if (utils::IntraComm::getSize() > 1) {
101 double surfaceArea = 0.0;
102 double surfaceAreaSum = 0.0;
103 utils::IntraComm::reduceSum(surfaceArea, surfaceAreaSum);
104 }
105 }
106}
107
108Eigen::VectorXd WatchIntegral::calculateIntegral(const mesh::PtrData &data) const
109{
110 int dim = data->getDimensions();
111 const Eigen::VectorXd &values = data->values();
112 Eigen::VectorXd sum = Eigen::VectorXd::Zero(dim);
113
114 if (_mesh->edges().empty() || (not _isScalingOn)) {
115 for (const auto &vertex : _mesh->vertices()) {
116 int offset = vertex.getID() * dim;
117 for (int i = 0; i < dim; i++) {
118 sum[i] += values[offset + i];
119 }
120 }
121 return sum;
122 } else { // Connectivity information is given
123 return mesh::integrateSurface(_mesh, data->values());
124 }
125}
126
128{
129 PRECICE_ASSERT(not _mesh->edges().empty());
130 double surfaceArea = 0.0;
131 if (_mesh->getDimensions() == 3) {
132 for (const auto &face : _mesh->triangles()) {
133 surfaceArea += face.getArea();
134 }
135 } else {
136 for (const auto &edge : _mesh->edges()) {
137 surfaceArea += edge.getLength();
138 }
139 }
140 return surfaceArea;
141}
142
143} // namespace precice::impl
#define PRECICE_WARN_IF(condition,...)
Definition LogMacros.hpp:21
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
void initialize()
Adds surface area information based on mesh connectivity.
std::vector< mesh::PtrData > _dataToExport
Eigen::VectorXd calculateIntegral(const mesh::PtrData &data) const
WatchIntegral(mesh::PtrMesh meshToWatch, const std::string &exportFilename, bool isScalingOn)
Constructor.
void exportIntegralData(double time)
Writes one line with data of the integral over the mesh into the output file.
io::TXTTableWriter _txtWriter
DataType
Constants defining possible data types to be written.
void addData(const std::string &name, DataType type)
Adds a data entry to the table.
void writeData(const std::string &name, int value)
Writes a integral scalar data value associated to the entry name.
static int getSize()
Number of ranks. This includes ranks from both participants, e.g. minimal size is 2.
Definition IntraComm.cpp:47
static bool isSecondary()
True if this process is running a secondary rank.
Definition IntraComm.cpp:57
static void reduceSum(precice::span< const double > sendData, precice::span< double > rcvData)
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
STL namespace.