Peano 4
Loading...
Searching...
No Matches
grid.cpp
Go to the documentation of this file.
1#include "grid.h"
2
3#include "GridVertex.h"
4#include "AutomatonState.h"
5#include "GridControlEvent.h"
6#include "GridStatistics.h"
7#include "Spacetree.h"
8#include "peano4/utils/Loop.h"
9
10namespace {
11 [[maybe_unused]] tarch::logging::Log _log("peano4::grid");
12} // namespace
13
14
15
17 switch (mode) {
19 return "aggressive-top-down";
20 break;
22 return "bottom-up";
23 break;
24 }
25 return "<undef>";
26}
27
29 std::ostringstream msg;
30 msg << "(" << toString(mode) << ",#cells=" << numberOfFineGridCells << ")";
31 return msg.str();
32}
33
34std::ostream& operator<<(std::ostream& out, const peano4::SplitInstruction& instruction) {
35 out << instruction.toString();
36 return out;
37}
38
39std::ostream& operator<<(std::ostream& out, const peano4::SplitInstruction::Mode& mode) {
41 return out;
42}
43
45 const peano4::grid::AutomatonState& state, const peano4::grid::GridControlEvent& event, double upscaleAutomatonState
46) {
47 assertion1(upscaleAutomatonState >= 1.0, upscaleAutomatonState);
48 return tarch::la::
49 allGreaterEquals(state.getX() + (0.5 - 0.5 * upscaleAutomatonState) * state.getH(), event.getOffset())
51 state.getX() + (0.5 + 0.5 * upscaleAutomatonState) * state.getH(), event.getOffset() + event.getWidth()
52 );
53}
54
59
64
66 [[maybe_unused]] GridVertex::State state,
67 [[maybe_unused]] const tarch::la::Vector<Dimensions,double>& x,
68 [[maybe_unused]] int level,
69 [[maybe_unused]] const tarch::la::Vector<TwoPowerD,int>& adjacentRanks,
70 [[maybe_unused]] bool isNewFineGridVertex
71) {
72 GridVertex result;
73
74 result.setState(state);
75 result.setAdjacentRanks(adjacentRanks);
76 result.setBackupOfAdjacentRanks(adjacentRanks);
77 result.setHasBeenAntecessorOfRefinedVertexInPreviousTreeSweep(not isNewFineGridVertex);
78 result.setIsAntecessorOfRefinedVertexInCurrentTreeSweep(not isNewFineGridVertex);
79#if PeanoDebug > 0
80 result.setX(x);
81#endif
82 result.setLevel(level);
83
84 // not required logically, but makes valgrind's memchecker happy
86
87 return result;
88}
89
90std::string peano4::grid::toString(const std::vector<GridControlEvent>& events) {
91 std::ostringstream msg;
92 msg << "(";
93 for (int i = 0; i < static_cast<int>(events.size()); i++) {
94 if (i > 0) {
95 msg << ",";
96 }
97 msg << events[i].toString();
98 }
99 msg << ")";
100 return msg.str();
101}
102
103std::string peano4::grid::toString(const std::list<GridControlEvent>& events) {
104 std::ostringstream msg;
105 msg << "[";
106 for (auto p : events) {
107 msg << p.toString();
108 }
109 msg << "]";
110 return msg.str();
111}
112
113std::vector<peano4::grid::GridControlEvent> peano4::grid::merge(
114 std::vector<GridControlEvent> events, [[maybe_unused]] const double Tolerance
115) {
116 logTraceInWith1Argument("merge(...)", events.size());
117
118 std::list<peano4::grid::GridControlEvent> refineEvents;
119 std::list<peano4::grid::GridControlEvent> eraseEvents;
120 for (auto p : events) {
121 if (p.getRefinementControl() == GridControlEvent::RefinementControl::Erase) {
122 eraseEvents.push_back(p);
123 } else {
124 refineEvents.push_back(p);
125 }
126 }
127 logDebug(
128 "merge(...)",
129 "have a total of "
130 << events.size() << " event(s) consisting of " << refineEvents.size() << " refine and " << eraseEvents.size()
131 << " erase events"
132 );
133
134 internal::removeEraseEventsThatAreCancelledByRefineEvents(refineEvents, eraseEvents);
135 logDebug(
136 "merge(...)",
137 "reduced to " << eraseEvents.size() << " erase event(s) by eliminating erase events that overlap with refinement"
138 );
139
140 internal::mergeAdjacentRefinementEvents(refineEvents, Tolerance);
141 internal::mergeAdjacentRefinementEvents(eraseEvents, Tolerance);
142 logDebug("merge(...)", "merged refine events into " << refineEvents.size() << " event(s)");
143
144 std::vector<peano4::grid::GridControlEvent> result;
145 result.insert(result.begin(), refineEvents.begin(), refineEvents.end());
146 result.insert(result.begin(), eraseEvents.begin(), eraseEvents.end());
147 logTraceOutWith1Argument("merge(...)", result.size());
148 return result;
149}
150
151void peano4::grid::internal::sort(std::list<peano4::grid::GridControlEvent>& events) {
152
153 auto compare = [](const peano4::grid::GridControlEvent& lhs, const peano4::grid::GridControlEvent& rhs) -> bool {
154 if ( tarch::la::smaller( tarch::la::volume(lhs.getH()), tarch::la::volume(rhs.getH()) ) ) return true;
155 for (int d = 0; d < Dimensions; d++) {
156 if (tarch::la::smaller(lhs.getOffset(d), rhs.getOffset(d)))
157 return true;
158 if (tarch::la::greater(lhs.getOffset(d), rhs.getOffset(d)))
159 return false;
160 }
161 return false;
162 };
163
164 events.sort(compare);
165}
166
168 const std::list<peano4::grid::GridControlEvent>& refineEvents, std::list<peano4::grid::GridControlEvent>& eraseEvents
169) {
170 std::list<peano4::grid::GridControlEvent>::iterator pEraseEvent = eraseEvents.begin();
171 while (pEraseEvent != eraseEvents.end()) {
172 bool overlaps = false;
173 for (auto refineEvent : refineEvents) {
174 overlaps |= refinementEventOverrulesCoarsening(refineEvent, *pEraseEvent);
175 }
176 if (overlaps) {
177 pEraseEvent = eraseEvents.erase(pEraseEvent);
178 } else {
179 pEraseEvent++;
180 }
181 }
182}
183
185 std::list<peano4::grid::GridControlEvent>& inputEvents, int Tolerance
186) {
187 bool hasMergedEvents = true;
188 while (hasMergedEvents) {
189 internal::sort(inputEvents);
190 hasMergedEvents = false;
191
192 std::list<peano4::grid::GridControlEvent>::iterator lhs = inputEvents.begin();
193 std::list<peano4::grid::GridControlEvent>::iterator rhs = inputEvents.begin();
194 rhs++;
195
196 while (rhs != inputEvents.end()) {
197 if (twoEventsAreAdjacent(*rhs, *lhs, Tolerance)) {
198 GridControlEvent newEvent = createBoundingBoxEvent(*rhs, *lhs);
199 logDebug(
200 "mergeAdjacentRefinementEvents(...)",
201 "merge two adjacent events "
202 << lhs->toString() << " and " << rhs->toString() << ") into " << newEvent.toString()
203 );
204 *lhs = newEvent;
205 lhs = inputEvents.erase(rhs);
206 rhs = lhs;
207 if (rhs != inputEvents.end())
208 rhs++;
209 hasMergedEvents = true;
210 } else {
211 lhs++;
212 rhs++;
213 }
214 }
215 }
216}
217
220 lhs.getNumberOfLocalUnrefinedCells() + rhs.getNumberOfLocalUnrefinedCells(),
221 lhs.getNumberOfRemoteUnrefinedCells() + rhs.getNumberOfRemoteUnrefinedCells(),
222 lhs.getNumberOfLocalRefinedCells() + rhs.getNumberOfLocalRefinedCells(),
223 lhs.getNumberOfRemoteRefinedCells() + rhs.getNumberOfRemoteRefinedCells(),
224 std::min(lhs.getStationarySweeps(), rhs.getStationarySweeps()),
225 lhs.getCoarseningHasBeenVetoed() or rhs.getCoarseningHasBeenVetoed(),
226 lhs.getRemovedEmptySubtree() or rhs.getRemovedEmptySubtree(),
227 tarch::la::min(lhs.getMinH(), rhs.getMinH())
228 );
229}
230
231void peano4::grid::clear(GridStatistics& statistics, bool isGlobalMasterTree) {
232 statistics.setNumberOfLocalUnrefinedCells(0);
234 statistics.setNumberOfLocalRefinedCells(isGlobalMasterTree ? 1 : 0);
235 statistics.setNumberOfRemoteRefinedCells(isGlobalMasterTree ? 0 : 1);
236 statistics.setCoarseningHasBeenVetoed(false);
237 statistics.setRemovedEmptySubtree(false);
238 statistics.setStationarySweeps(statistics.getStationarySweeps() + 1);
239 statistics.setMinH(tarch::la::Vector<Dimensions, double>(std::numeric_limits<double>::max()));
240}
241
243 bool result = false;
244 dfor2(k) result |= willVertexBeRefined(vertices[kScalar]);
245 result |= hasVertexBeenRefined(vertices[kScalar]);
246 enddforx return result;
247}
248
250 return vertex.getState() == GridVertex::State::Refining or vertex.getState() == GridVertex::State::Refined
251 or vertex.getState() == GridVertex::State::EraseTriggered;
252}
253
255 return vertex.getState() == GridVertex::State::Refined or vertex.getState() == GridVertex::State::EraseTriggered
256 or vertex.getState() == GridVertex::State::Erasing;
257}
258
260 std::bitset<TwoPowerD> bitset;
261 for (int i = 0; i < TwoPowerD; i++) {
262 assertion(not willVertexBeRefined(vertices[i]) or vertices[i].getState() != GridVertex::State::HangingVertex);
263 bitset.set(i, willVertexBeRefined(vertices[i]));
264 }
265 return bitset;
266}
267
269 std::bitset<TwoPowerD> bitset;
270 for (int i = 0; i < TwoPowerD; i++) {
271 assertion(not hasVertexBeenRefined(vertices[i]) or vertices[i].getState() != GridVertex::State::HangingVertex);
272 bitset.set(i, hasVertexBeenRefined(vertices[i]));
273 }
274 return bitset;
275}
276
278 switch (type) {
279 case VertexType::New:
280 return "new";
281 case VertexType::Hanging:
282 return "hanging";
283 case VertexType::Persistent:
284 return "persistent";
285 case VertexType::Delete:
286 return "delete";
287 }
288 return "<undef>";
289}
290
292 switch (type) {
293 case FaceType::New:
294 return "new";
295 case FaceType::Hanging:
296 return "hanging";
297 case FaceType::Persistent:
298 return "persistent";
299 case FaceType::Delete:
300 return "delete";
301 }
302 return "<undef>";
303}
304
306 switch (type) {
307 case CellType::New:
308 return "new";
309 case CellType::Persistent:
310 return "persistent";
311 case CellType::Delete:
312 return "delete";
313 }
314 return "<undef>";
315}
316
318 switch (state) {
319 case SpacetreeState::EmptyRun:
320 return "empty-run";
321 case SpacetreeState::NewRoot:
322 return "new-root";
323 case SpacetreeState::NewFromSplit:
324 return "new-from-split";
325 case SpacetreeState::Running:
326 return "running";
327 case SpacetreeState::JoinTriggered:
328 return "join-triggered";
329 case SpacetreeState::Joining:
330 return "joining";
331 case SpacetreeState::Joined:
332 return "joined";
333 }
334 return "<undef>";
335}
336
339) {
340 bool twoEventsOverlap = true;
341 for (int d = 0; d < Dimensions; d++) {
342 twoEventsOverlap &= lhs.getOffset()(d) + lhs.getWidth()(d) >= rhs.getOffset()(d);
343 twoEventsOverlap &= lhs.getOffset()(d) <= rhs.getOffset()(d) + rhs.getWidth()(d);
344 }
345 return twoEventsOverlap;
346};
347
350) {
351 bool sameType
352 = (lhs.getRefinementControl() == GridControlEvent::RefinementControl::Erase
353 and rhs.getRefinementControl() == GridControlEvent::RefinementControl::Erase)
354 or (lhs.getRefinementControl() == GridControlEvent::RefinementControl::Refine and rhs.getRefinementControl() == GridControlEvent::RefinementControl::Refine);
355
356 return sameType and tarch::la::equals(lhs.getOffset(), rhs.getOffset())
357 and tarch::la::equals(lhs.getWidth(), rhs.getWidth()) and tarch::la::equals(lhs.getH(), rhs.getH());
358};
359
361 const peano4::grid::GridControlEvent& refineEvent, const peano4::grid::GridControlEvent& eraseEvent
362) {
363 return twoEventsOverlap(refineEvent, eraseEvent)
364 and refineEvent.getRefinementControl() == GridControlEvent::RefinementControl::Refine
365 and eraseEvent.getRefinementControl() == GridControlEvent::RefinementControl::Erase
366 and tarch::la::oneSmallerEquals(1.0 / 3.0 * refineEvent.getH(), eraseEvent.getH());
367};
368
370 const peano4::grid::GridControlEvent& lhs, const peano4::grid::GridControlEvent& rhs, double tolerance
371) {
372 tarch::la::Vector<Dimensions, double> boundingEventOffset = tarch::la::min(lhs.getOffset(), rhs.getOffset());
374 boundingEventSize = tarch::la::max(lhs.getOffset() + lhs.getWidth(), rhs.getOffset() + rhs.getWidth())
375 - boundingEventOffset;
376
377 const double relativeTolerance = tarch::la::relativeEps(
378 tarch::la::volume(lhs.getWidth()), tarch::la::volume(rhs.getWidth()), tolerance
379 );
380 const double fusedEventsSize = tarch::la::volume(lhs.getWidth()) + tarch::la::volume(rhs.getWidth());
381 return tarch::la::equals(lhs.getH(), rhs.getH(), relativeTolerance)
382 and tarch::la::smallerEquals(tarch::la::volume(boundingEventSize), fusedEventsSize, relativeTolerance);
383};
384
387) {
388 tarch::la::Vector<Dimensions, double> boundingEventOffset = tarch::la::min(lhs.getOffset(), rhs.getOffset());
390 boundingEventSize = tarch::la::max(lhs.getOffset() + lhs.getWidth(), rhs.getOffset() + rhs.getWidth())
391 - boundingEventOffset;
392 return GridControlEvent(rhs.getRefinementControl(), boundingEventOffset, boundingEventSize, rhs.getH());
393};
394
395void peano4::grid::reduceGridControlEvents([[maybe_unused]] std::vector<GridControlEvent>& events) {
396#ifdef Parallel
397 std::vector<GridControlEvent> localEvents = events;
398 int localSize = events.size();
399 int* sizes = new int[tarch::mpi::Rank::getInstance().getNumberOfRanks()];
400
401 logDebug("reduceGridControlEvents(...)", "feed local data of size " << localEvents.size() << " into allgather");
402
403 MPI_Allgather(&localSize, 1, MPI_INT, sizes, 1, MPI_INT, tarch::mpi::Rank::getInstance().getCommunicator());
404
405 int totalSize = 0;
406 for (int i = 0; i < tarch::mpi::Rank::getInstance().getNumberOfRanks(); i++) {
407 totalSize += sizes[i];
408 }
409 events.resize(totalSize);
410 logDebug(
411 "reduceGridControlEvents(...)",
412 "we have "
413 << totalSize << " grid control event(s) over all ranks, while the local rank hosts " << events.size()
414 << " event(s)"
415 );
416
417 int offset = 0;
418 for (int i = 0; i < tarch::mpi::Rank::getInstance().getNumberOfRanks(); i++) {
420 for (int j = 0; j < sizes[i]; j++) {
421 events[offset + j] = localEvents[j];
422 }
423 }
424
425 if (sizes[i] > 0) {
426 MPI_Bcast(
427 events.data() + offset,
428 sizes[i],
429 GridControlEvent::getGlobalCommunciationDatatype(),
430 i,
432 );
433 offset += sizes[i];
434 }
435 }
436
437 delete[] sizes;
438#endif
439}
#define assertion1(expr, param)
#define assertion(expr)
std::ostream & operator<<(std::ostream &out, const peano4::datamanagement::CellMarker &marker)
Definition CellMarker.cpp:7
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
#define TwoPowerD
Definition Globals.h:19
AutomatonState state
#define logDebug(methodName, logMacroMessageStream)
Definition Log.h:50
#define logTraceOutWith1Argument(methodName, argument0)
Definition Log.h:380
#define logTraceInWith1Argument(methodName, argument0)
Definition Log.h:370
#define dfor2(counter)
Shortcut For dfor(counter,2)
Definition Loop.h:906
#define enddforx
I prefer to use this macro for dforx instead of a closing bracket as many syntax parser fail otherwis...
Definition Loop.h:928
Log Device.
Definition Log.h:516
int getNumberOfRanks() const
Definition Rank.cpp:551
static Rank & getInstance()
This operation returns the singleton instance.
Definition Rank.cpp:538
int getRank() const
Return rank of this node.
Definition Rank.cpp:528
MPI_Comm getCommunicator() const
Definition Rank.cpp:544
std::string toString(Filter filter)
Definition convert.cpp:170
peano4::grid::GridStatistics operator+(peano4::grid::GridStatistics lhs, peano4::grid::GridStatistics rhs)
Definition grid.cpp:218
bool twoEventsAreAdjacent(const peano4::grid::GridControlEvent &lhs, const peano4::grid::GridControlEvent &rhs, double Tolerance)
Are two events adjacent.
Definition grid.cpp:369
void mergeAdjacentRefinementEvents(std::list< peano4::grid::GridControlEvent > &inputEvents, int Tolerance)
Merge adjacent events.
Definition grid.cpp:184
void removeEraseEventsThatAreCancelledByRefineEvents(const std::list< peano4::grid::GridControlEvent > &refineEvents, std::list< peano4::grid::GridControlEvent > &eraseEvents)
This is the first thing I do.
Definition grid.cpp:167
bool refinementEventOverrulesCoarsening(const peano4::grid::GridControlEvent &refineEvent, const peano4::grid::GridControlEvent &eraseEvent)
A refinement event overrules the coarsening if.
Definition grid.cpp:360
peano4::grid::GridControlEvent createBoundingBoxEvent(const peano4::grid::GridControlEvent &lhs, const peano4::grid::GridControlEvent &rhs)
Definition grid.cpp:385
void sort(std::list< peano4::grid::GridControlEvent > &events)
Sort grid control events geometrically.
Definition grid.cpp:151
bool twoEventsOverlap(const peano4::grid::GridControlEvent &lhs, const peano4::grid::GridControlEvent &rhs)
Helper function which helps us throughout the merge.
Definition grid.cpp:337
bool equals(const peano4::grid::GridControlEvent &lhs, const peano4::grid::GridControlEvent &rhs)
Definition grid.cpp:348
bool isSpacetreeNodeRefined(GridVertex vertices[TwoPowerD])
A spacetree node is refined if any of its adjacent vertices holds one of the following flags:
Definition grid.cpp:242
std::string toString(VertexType type)
Definition grid.cpp:277
void clear(GridStatistics &statistics, bool isGlobalMasterTree)
The term clear() is not 100% correct, as the number of stationary traversals is not reset to a dummy ...
Definition grid.cpp:231
void reduceGridControlEvents(std::vector< GridControlEvent > &events)
Peano 4 does not reduce any grid control events globally.
Definition grid.cpp:395
GridVertex createVertex(GridVertex::State state, const tarch::la::Vector< Dimensions, double > &x, int level, const tarch::la::Vector< TwoPowerD, int > &adjacentRanks, bool isNewFineGridVertex)
Factory mechanism.
Definition grid.cpp:65
std::bitset< TwoPowerD > haveVerticesBeenRefined(GridVertex vertices[TwoPowerD])
Definition grid.cpp:268
bool overlaps(const tarch::la::Vector< Dimensions, double > &x, const GridControlEvent &event)
Definition grid.cpp:60
std::bitset< TwoPowerD > willVerticesBeRefined(GridVertex vertices[TwoPowerD])
A vertex is unrefined if it is hanging.
Definition grid.cpp:259
bool hasVertexBeenRefined(const GridVertex &vertex)
A vertex has been refined if it is (already) refined or is erasing or the erase has been triggered.
Definition grid.cpp:254
std::vector< GridControlEvent > merge(std::vector< GridControlEvent > events, const double Tolerance=0.1)
Merge set of refinement/coarsening commands.
Definition grid.cpp:113
bool willVertexBeRefined(const GridVertex &vertex)
A vertex will be refined if it is already refined or currently refining.
Definition grid.cpp:249
bool isContained(const AutomatonState &x, const GridControlEvent &event, double upscaleAutomatonState=1.0)
isContained() is defined over the closed interval, i.e.
Definition grid.cpp:44
Scalar volume(const Vector< Size, Scalar > &vector)
Computes the volume of the tetrahedron spanned by the Cartesian unit vectors scaled by the correspond...
bool allGreaterEquals(const Vector< Size, Scalar > &lhs, const Scalar &cmp, const Scalar tolerance=NUMERICAL_ZERO_DIFFERENCE)
static bool smaller(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE) InlineMethod
Smaller operator for floating point values.
bool greater(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
bool oneSmallerEquals(const Vector< Size, Scalar > &lhs, const Scalar &cmp, const Scalar tolerance=NUMERICAL_ZERO_DIFFERENCE)
double relativeEps(double valueA, double valueB=std::numeric_limits< double >::min(), double eps=NUMERICAL_ZERO_DIFFERENCE)
Determine a relative tolerance from one or two values.
Definition Scalar.cpp:14
double max(double a, double b, double c)
I need the maximum of three values all the time, to I decided to write a function for this.
Definition Scalar.cpp:8
bool smallerEquals(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
bool allSmallerEquals(const Vector< Size, Scalar > &lhs, const Scalar &cmp, const Scalar tolerance=NUMERICAL_ZERO_DIFFERENCE)
bool equals(const Matrix< Rows, Cols, Scalar > &lhs, const Matrix< Rows, Cols, Scalar > &rhs, const Scalar &tolerance=NUMERICAL_ZERO_DIFFERENCE)
Compares to matrices on equality by means of a numerical accuracy.
Scalar min(const Vector< Size, Scalar > &vector)
Returns the element with minimal value (NOT absolute value).
tarch::logging::Log _log("examples::unittests")
Instruction to split.
Definition grid.h:34
static std::string toString(Mode mode)
Definition grid.cpp:16
std::string toString() const
Definition grid.cpp:28
tarch::la::Vector< Dimensions, double > getWidth() const
peano4::grid::GridControlEvent::RefinementControl getRefinementControl() const
tarch::la::Vector< Dimensions, double > getOffset() const
tarch::la::Vector< Dimensions, double > getH() const
void setNumberOfLocalRefinedCells(int value)
void setNumberOfRemoteUnrefinedCells(int value)
void setMinH(const tarch::la::Vector< Dimensions, double > &value)
void setCoarseningHasBeenVetoed(bool value)
void setNumberOfLocalUnrefinedCells(int value)
void setNumberOfRemoteRefinedCells(int value)
tarch::la::Vector< Dimensions, double > getMinH() const
void setRemovedEmptySubtree(bool value)
void setBackupOfAdjacentRanks(const tarch::la::Vector< TwoPowerD, int > &value)
void setIsAntecessorOfRefinedVertexInCurrentTreeSweep(bool value)
void setLevel(int value)
void setNumberOfAdjacentRefinedLocalCells(int value)
void setHasBeenAntecessorOfRefinedVertexInPreviousTreeSweep(bool value)
void setState(State value)
void setAdjacentRanks(const tarch::la::Vector< TwoPowerD, int > &value)
peano4::grid::GridVertex::State getState() const
Simple vector class.
Definition Vector.h:134