1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
|
#include "WarpX.H"
#include <AMReX_Config.H>
#ifdef AMREX_USE_EB
# include <AMReX_EB2.H>
# include <AMReX_ParmParse.H>
#endif
void
WarpX::InitEB ()
{
#ifdef AMREX_USE_EB
BL_PROFILE("InitEB");
amrex::ParmParse pp_eb2("eb2");
if (!pp_eb2.contains("geom_type")) {
std::string geom_type = "all_regular";
pp_eb2.add("geom_type", geom_type); // use all_regular by default
}
amrex::EB2::Build(Geom(maxLevel()), maxLevel(), maxLevel());
#endif
}
/**
* \brief Compute the length of the mesh edges. Here the length is a value in [0, 1].
* An edge of length 0 is fully covered.
*/
void
WarpX::ComputeEdgeLengths () {
#ifdef AMREX_USE_EB
BL_PROFILE("ComputeEdgeLengths");
auto const eb_fact = fieldEBFactory(maxLevel());
auto const &flags = eb_fact.getMultiEBCellFlagFab();
auto const &edge_centroid = eb_fact.getEdgeCent();
for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi){
amrex::Box const &box = mfi.validbox();
amrex::FabType fab_type = flags[mfi].getType(box);
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim){
auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi);
if (fab_type == amrex::FabType::regular) {
// every cell in box is all regular
amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge_lengths_dim).ixType()),
[=](int i, int j, int k) {
edge_lengths_dim(i, j, k) = 1.;
});
} else if (fab_type == amrex::FabType::covered) {
// every cell in box is all covered
amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge_lengths_dim).ixType()),
[=](int i, int j, int k) {
edge_lengths_dim(i, j, k) = 0.;
});
} else {
auto const &edge_cent = edge_centroid[idim]->const_array(mfi);
amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge_cent).ixType()),
[=](int i, int j, int k) {
if (edge_cent(i, j, k) == amrex::Real(-1.0)) {
// This edge is all covered
edge_lengths_dim(i, j, k) = 0.;
} else if (edge_cent(i, j, k) == amrex::Real(1.0)) {
// This edge is all open
edge_lengths_dim(i, j, k) = 1.;
} else {
// This edge is cut.
edge_lengths_dim(i, j, k) = 1 - abs(amrex::Real(2.0)*edge_cent(i, j, k));
}
});
}
}
}
#endif
}
/**
* \brief Compute the ara of the mesh faces. Here the area is a value in [0, 1].
* An edge of area 0 is fully covered.
*/
void
WarpX::ComputeFaceAreas () {
#ifdef AMREX_USE_EB
BL_PROFILE("ComputeFaceAreas");
auto const eb_fact = fieldEBFactory(maxLevel());
auto const &flags = eb_fact.getMultiEBCellFlagFab();
auto const &area_frac = eb_fact.getAreaFrac();
for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) {
amrex::Box const &box = mfi.validbox();
amrex::FabType fab_type = flags[mfi].getType(box);
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi);
if (fab_type == amrex::FabType::regular) {
// every cell in box is all regular
amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face_areas_dim).ixType()),
[=](int i, int j, int k) {
face_areas_dim(i, j, k) = amrex::Real(1.);
});
} else if (fab_type == amrex::FabType::covered) {
// every cell in box is all covered
amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face_areas_dim).ixType()),
[=](int i, int j, int k) {
face_areas_dim(i, j, k) = amrex::Real(0.);
});
} else {
auto const &face = area_frac[idim]->const_array(mfi);
amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face).ixType()),
[=](int i, int j, int k) {
face_areas_dim(i, j, k) = face(i, j, k);
});
}
}
}
#endif
}
/**
* \brief Scale the edges lengths by the mesh width to obtain the real lengths.
*/
void
WarpX::ScaleEdges () {
#ifdef AMREX_USE_EB
BL_PROFILE("ScaleEdges");
auto const &cell_size = CellSize(maxLevel());
auto const eb_fact = fieldEBFactory(maxLevel());
auto const &flags = eb_fact.getMultiEBCellFlagFab();
for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) {
amrex::Box const &box = mfi.validbox();
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi);
amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge_lengths_dim).ixType()),
[=](int i, int j, int k) {
edge_lengths_dim(i, j, k) *= cell_size[idim];
});
}
}
#endif
}
/**
* \brief Scale the edges areas by the mesh width to obtain the real areas.
*/
void
WarpX::ScaleAreas() {
#ifdef AMREX_USE_EB
BL_PROFILE("ScaleAreas");
auto const& cell_size = CellSize(maxLevel());
amrex::Real full_area;
auto const eb_fact = fieldEBFactory(maxLevel());
auto const &flags = eb_fact.getMultiEBCellFlagFab();
for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) {
amrex::Box const &box = mfi.validbox();
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (idim == 0) {
full_area = cell_size[1]*cell_size[2];
} else if (idim == 1) {
full_area = cell_size[0]*cell_size[2];
} else {
full_area = cell_size[0]*cell_size[1];
}
auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi);
amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face_areas_dim).ixType()),
[=](int i, int j, int k) {
face_areas_dim(i, j, k) *= full_area;
});
}
}
#endif
}
|
