Skip to content

Commit e5616c1

Browse files
cwsmithcwsmith
committed
formatting
1 parent 9026150 commit e5616c1

File tree

3 files changed

+206
-204
lines changed

3 files changed

+206
-204
lines changed

‎src/MeshField.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -313,7 +313,7 @@ class OmegahMeshField {
313313
auto field_view = field.vtxField.serialize();
314314
Omega_h::Write<FieldDataType> field_write(field_view);
315315
mesh.add_tag(0, "field", 1, Omega_h::read(field_write), false,
316-
Omega_h::ArrayType::VectorND);
316+
Omega_h::ArrayType::VectorND);
317317
Omega_h::vtk::write_parallel("foo.vtk", &mesh, mesh.dim());
318318
}
319319

‎src/MeshField_SPR_ErrorEstimator.hpp‎

Lines changed: 117 additions & 127 deletions
Original file line numberDiff line numberDiff line change
@@ -1,52 +1,45 @@
11
#ifndef THWAITES_ERROR_ESTIMATOR_H
22
#define THWAITES_ERROR_ESTIMATOR_H
3-
#include <Omega_h_mesh.hpp>
43
#include <KokkosController.hpp>
54
#include <MeshField.hpp>
65
#include <MeshField_Element.hpp>
7-
#include <MeshField_Integrate.hpp>
86
#include <MeshField_Fail.hpp>
97
#include <MeshField_For.hpp>
8+
#include <MeshField_Integrate.hpp>
109
#include <MeshField_ShapeField.hpp>
11-
10+
#include <Omega_h_mesh.hpp>
1211

1312
/* useful for initializing values to quickly
1413
detect "uninitialized" value bugs */
15-
static double getNaN()
16-
{
17-
return std::numeric_limits<double>::quiet_NaN();
18-
}
14+
static double getNaN() { return std::numeric_limits<double>::quiet_NaN(); }
1915

2016
/* common base for Scalar Integrator. */
2117
class SInt : public MeshField::Integrator {
22-
public:
23-
SInt(int order):
24-
MeshField::Integrator(order),r(0)
25-
{}
26-
void reset() {r=0;}
27-
MeshField::Real r;
18+
public:
19+
SInt(int order) : MeshField::Integrator(order), r(0) {}
20+
void reset() { r = 0; }
21+
MeshField::Real r;
2822
};
2923

30-
template <typename ShapeField>
31-
class Estimation {
32-
private:
24+
template <typename ShapeField> class Estimation {
25+
private:
3326
Estimation() {}
34-
public:
35-
Estimation(Omega_h::Mesh& mesh_in, Omega_h::Reals eps_in, ShapeField& eps_star_in,
36-
MeshField::Real tolerance_in)
37-
: mesh(mesh_in), eps(eps_in), eps_star(eps_star_in),
38-
tolerance(tolerance_in), size_factor(getNaN())
39-
{
27+
28+
public:
29+
Estimation(Omega_h::Mesh &mesh_in, Omega_h::Reals eps_in,
30+
ShapeField &eps_star_in, MeshField::Real tolerance_in)
31+
: mesh(mesh_in), eps(eps_in), eps_star(eps_star_in),
32+
tolerance(tolerance_in), size_factor(getNaN()) {
4033
/* note that getOrder being used to convey this meaning
4134
is a bit of a hack, but looks decent if you don't
4235
try to define the FieldShape API too rigorously */
43-
integration_order = 1; //FIXME get from eps_star
36+
integration_order = 1; // FIXME get from eps_star
4437
/* so far recovery order is directly tied to the
4538
mesh's coordinate field order, coordinate this
4639
with field recovery code */
47-
recovered_order = 1; //FIXME
40+
recovered_order = 1; // FIXME
4841
}
49-
Omega_h::Mesh& mesh;
42+
Omega_h::Mesh &mesh;
5043
/* the maximum polynomial order that can be
5144
integrated exactly using the input integration points.
5245
not necessarily equal to recovered_order, sometimes
@@ -61,7 +54,7 @@ class Estimation {
6154
/* the recovered field, consisting of nodal values
6255
of the quantity of interest.
6356
this uses a Lagrange basis of recovered_order */
64-
ShapeField& eps_star;
57+
ShapeField &eps_star;
6558
/* the acceptable margin of error, expressed as a factor
6659
greater than zero.
6760
setting this equal to zero would request zero element
@@ -76,168 +69,165 @@ class Estimation {
7669
desired element size = current size * current error * size_factor */
7770
Omega_h::Real size_factor;
7871
/* a temporary field storing desired sizes at elements */
79-
Kokkos::View<MeshField::Real*> element_size;
72+
Kokkos::View<MeshField::Real *> element_size;
8073
/* the resulting size field, recovered from the element_size field
8174
(using a local average recovery method much weaker than SPR) */
8275
Omega_h::Reals size;
8376
};
8477

8578
/* computes $\|f\|^2$ */
86-
template<typename EstimationT, typename OmegahMeshField>
87-
class SelfProduct : public SInt
88-
{
89-
public:
90-
SelfProduct(EstimationT& estimation_in, OmegahMeshField& omf_in):
91-
SInt(estimation_in.integration_order),
92-
estimation(estimation_in),
93-
omf(omf_in)
94-
{
95-
}
96-
void atPoints(Kokkos::View<MeshField::Real**> p,
97-
Kokkos::View<MeshField::Real*> w,
98-
Kokkos::View<MeshField::Real*> dV)
99-
{
100-
101-
std::cerr << "SelfProduct::atPoints(...)\n";
102-
const size_t numPtsPerElem = p.extent(0)/estimation.mesh.nelems();
103-
auto eps_star_atPts = omf.triangleLocalPointEval(p,numPtsPerElem,
104-
estimation.eps_star);
105-
106-
r = 0;
107-
Kokkos::parallel_reduce(
79+
template <typename EstimationT, typename OmegahMeshField>
80+
class SelfProduct : public SInt {
81+
public:
82+
SelfProduct(EstimationT &estimation_in, OmegahMeshField &omf_in)
83+
: SInt(estimation_in.integration_order), estimation(estimation_in),
84+
omf(omf_in) {}
85+
void atPoints(Kokkos::View<MeshField::Real **> p,
86+
Kokkos::View<MeshField::Real *> w,
87+
Kokkos::View<MeshField::Real *> dV) {
88+
89+
std::cerr << "SelfProduct::atPoints(...)\n";
90+
const size_t numPtsPerElem = p.extent(0) / estimation.mesh.nelems();
91+
auto eps_star_atPts =
92+
omf.triangleLocalPointEval(p, numPtsPerElem, estimation.eps_star);
93+
94+
r = 0;
95+
Kokkos::parallel_reduce(
10896
"eval", estimation.mesh.nelems(),
10997
KOKKOS_LAMBDA(const int &elm, MeshField::Real &r_local) {
11098
const auto first = elm * numPtsPerElem;
11199
const auto last = first + numPtsPerElem;
112100
for (auto pt = first; pt < last; pt++) {
113-
const auto vPt = eps_star_atPts(pt,0);
101+
const auto vPt = eps_star_atPts(pt, 0);
114102
const auto wPt = w(pt);
115103
const auto dVPt = dV(pt);
116104
r_local += (vPt * vPt) * wPt * dVPt;
117105
}
118106
},
119107
r);
108+
}
120109

121-
}
122-
private:
123-
EstimationT& estimation;
124-
OmegahMeshField& omf;
110+
private:
111+
EstimationT &estimation;
112+
OmegahMeshField &omf;
125113
};
126114

127115
/* computes the integral over the element of the
128116
sum of the squared differences between the
129117
original and recovered fields */
130-
template<typename EstimationT, typename OmegahMeshField>
131-
class Error : public SInt
132-
{
133-
public:
134-
Error(EstimationT& estimation_in, OmegahMeshField& omf_in):
135-
SInt(estimation_in.integration_order),
136-
estimation(estimation_in),
137-
omf(omf_in),
138-
errorNorm("errorNorm", estimation_in.mesh.nelems())
139-
{
140-
}
141-
void atPoints(Kokkos::View<MeshField::Real**> p,
142-
Kokkos::View<MeshField::Real*> w,
143-
Kokkos::View<MeshField::Real*> dV)
144-
{
145-
std::cerr << "SelfProduct::atPoints(...)\n";
146-
const size_t numPtsPerElem = p.extent(0)/estimation.mesh.nelems();
147-
//FIXME eps isn't a ShapeField so we can't call
148-
// omf.triangleLocalPointEval. For now, just get
149-
// the value from the element and assert that there is
150-
// one integration point per element.
151-
assert(numPtsPerElem == 1);
152-
auto eps_star_atPts = omf.triangleLocalPointEval(p,numPtsPerElem,
153-
estimation.eps_star);
154-
double meshDim = estimation.mesh.dim();
155-
double orderP = estimation.recovered_order;
156-
157-
const auto& eps = estimation.eps;
158-
r = 0;
159-
Kokkos::parallel_reduce(
118+
template <typename EstimationT, typename OmegahMeshField>
119+
class Error : public SInt {
120+
public:
121+
Error(EstimationT &estimation_in, OmegahMeshField &omf_in)
122+
: SInt(estimation_in.integration_order), estimation(estimation_in),
123+
omf(omf_in), errorNorm("errorNorm", estimation_in.mesh.nelems()) {}
124+
void atPoints(Kokkos::View<MeshField::Real **> p,
125+
Kokkos::View<MeshField::Real *> w,
126+
Kokkos::View<MeshField::Real *> dV) {
127+
std::cerr << "SelfProduct::atPoints(...)\n";
128+
const size_t numPtsPerElem = p.extent(0) / estimation.mesh.nelems();
129+
// FIXME eps isn't a ShapeField so we can't call
130+
// omf.triangleLocalPointEval. For now, just get
131+
// the value from the element and assert that there is
132+
// one integration point per element.
133+
assert(numPtsPerElem == 1);
134+
auto eps_star_atPts =
135+
omf.triangleLocalPointEval(p, numPtsPerElem, estimation.eps_star);
136+
double meshDim = estimation.mesh.dim();
137+
double orderP = estimation.recovered_order;
138+
139+
const auto &eps = estimation.eps;
140+
r = 0;
141+
Kokkos::parallel_reduce(
160142
"eval", estimation.mesh.nelems(),
161143
KOKKOS_CLASS_LAMBDA(const int &elm, MeshField::Real &r_local) {
162144
const auto first = elm * numPtsPerElem;
163145
const auto last = first + numPtsPerElem;
164146
const auto eps_elm = eps[elm];
165147
MeshField::Real sum = 0;
166148
for (auto pt = first; pt < last; pt++) {
167-
const auto eps_star_Pt = eps_star_atPts(pt,0);
149+
const auto eps_star_Pt = eps_star_atPts(pt, 0);
168150
const auto diff = eps_elm - eps_star_Pt;
169151
const auto wPt = w(pt);
170152
const auto dVPt = dV(pt);
171153
sum += (diff * diff) * wPt * dVPt;
172154
}
173-
assert(sum>0);
174-
r_local += Kokkos::pow(Kokkos::sqrt(sum), ((2 * meshDim) / (2 * orderP + meshDim)));
175-
errorNorm(elm) = Kokkos::pow(Kokkos::sqrt(sum), -(2 / (2 * orderP + meshDim)));
155+
assert(sum > 0);
156+
r_local += Kokkos::pow(Kokkos::sqrt(sum),
157+
((2 * meshDim) / (2 * orderP + meshDim)));
158+
errorNorm(elm) =
159+
Kokkos::pow(Kokkos::sqrt(sum), -(2 / (2 * orderP + meshDim)));
176160
},
177-
r); // $\sum_{i=1}^n \|e_\epsilon\|^{\frac{2d}{2p+d}}$
178-
}
179-
EstimationT& estimation;
180-
OmegahMeshField& omf;
181-
Kokkos::View<MeshField::Real*> errorNorm; // (||e_eps||_e)^(-2/(2p+d))
161+
r); // $\sum_{i=1}^n \|e_\epsilon\|^{\frac{2d}{2p+d}}$
162+
}
163+
EstimationT &estimation;
164+
OmegahMeshField &omf;
165+
Kokkos::View<MeshField::Real *> errorNorm; // (||e_eps||_e)^(-2/(2p+d))
182166
};
183167

184-
//TODO move this into Estimation class
185-
template<typename EstimationT, typename OmegahMeshField, typename FieldElement, typename ErrorT>
186-
void computeSizeFactor(EstimationT& e, OmegahMeshField& omf, FieldElement& coordFe, ErrorT& errorIntegrator) {
187-
SelfProduct sp(e,omf);
168+
// TODO move this into Estimation class
169+
template <typename EstimationT, typename OmegahMeshField, typename FieldElement,
170+
typename ErrorT>
171+
void computeSizeFactor(EstimationT &e, OmegahMeshField &omf,
172+
FieldElement &coordFe, ErrorT &errorIntegrator) {
173+
SelfProduct sp(e, omf);
188174
sp.process(coordFe);
189175
const double epsStarNorm = Kokkos::sqrt(sp.r);
190176
std::cout << "SelfProduct: " << epsStarNorm << "\n";
191177
const double a = e.tolerance * e.tolerance * // (n^hat)^2
192-
epsStarNorm * epsStarNorm; // ||e*||^2
193-
const double b = a / errorIntegrator.r; // term in parenthesis in section 4 of spr.tex
178+
epsStarNorm * epsStarNorm; // ||e*||^2
179+
const double b =
180+
a / errorIntegrator.r; // term in parenthesis in section 4 of spr.tex
194181
const double p = e.recovered_order;
195182
e.size_factor = Kokkos::pow(b, 1.0 / (2.0 * p));
196183
std::cout << "size_factor: " << e.size_factor << "\n";
197-
198184
}
199185

200186
/* computes h_e^new from section 4 of spr.tex */
201-
template<typename EstimationT, typename ErrorT>
202-
void getElementSizeField(EstimationT& e, ErrorT& errorIntegrator) {
203-
Kokkos::View<MeshField::Real*> eSize("eSize", e.mesh.nelems());
187+
template <typename EstimationT, typename ErrorT>
188+
void getElementSizeField(EstimationT &e, ErrorT &errorIntegrator) {
189+
Kokkos::View<MeshField::Real *> eSize("eSize", e.mesh.nelems());
204190
const auto errorNorm = errorIntegrator.errorNorm;
205191
const auto size_factor = e.size_factor;
206192
const auto currentElmSize = e.mesh.ask_sizes();
207-
e.mesh.template add_tag<MeshField::Real>(e.mesh.dim(), "curElmSize", 1, currentElmSize,
208-
false, Omega_h::ArrayType::VectorND);
209-
Kokkos::parallel_for(e.mesh.nelems(), KOKKOS_LAMBDA(const int elm) {
210-
const double h = currentElmSize[elm]; // h_e^current
211-
eSize(elm) = h * errorNorm(elm) * size_factor;
212-
});
193+
e.mesh.template add_tag<MeshField::Real>(e.mesh.dim(), "curElmSize", 1,
194+
currentElmSize, false,
195+
Omega_h::ArrayType::VectorND);
196+
Kokkos::parallel_for(
197+
e.mesh.nelems(), KOKKOS_LAMBDA(const int elm) {
198+
const double h = currentElmSize[elm]; // h_e^current
199+
eSize(elm) = h * errorNorm(elm) * size_factor;
200+
});
213201
e.element_size = eSize;
214202
}
215203

216-
Kokkos::View<MeshField::Real*>
217-
averageToVertex(Omega_h::Mesh& mesh, const Kokkos::View<MeshField::Real*>& elmSize) {
204+
Kokkos::View<MeshField::Real *>
205+
averageToVertex(Omega_h::Mesh &mesh,
206+
const Kokkos::View<MeshField::Real *> &elmSize) {
218207
Omega_h::Write<MeshField::Real> elmSize_oh(elmSize);
219-
mesh.add_tag<MeshField::Real>(mesh.dim(), "elmSizeField", 1, elmSize_oh, false,
220-
Omega_h::ArrayType::VectorND);
221-
Kokkos::View<MeshField::Real*> sizeField("sizeField", mesh.nverts());
208+
mesh.add_tag<MeshField::Real>(mesh.dim(), "elmSizeField", 1, elmSize_oh,
209+
false, Omega_h::ArrayType::VectorND);
210+
Kokkos::View<MeshField::Real *> sizeField("sizeField", mesh.nverts());
222211
const auto v2e = mesh.ask_up(Omega_h::VERT, mesh.dim());
223212
const auto v2e_offsets = v2e.a2ab;
224213
const auto v2e_values = v2e.ab2b;
225-
Kokkos::parallel_for(mesh.nverts(), KOKKOS_LAMBDA(const int vtx) {
226-
MeshField::Real s = 0;
227-
for (auto idx = v2e_offsets[vtx]; idx < v2e_offsets[vtx + 1]; ++idx) {
228-
const auto elm = v2e_values[idx];
229-
s += elmSize(elm);
230-
}
231-
const auto numUpElms = v2e_offsets[vtx+1] - v2e_offsets[vtx];
232-
sizeField(vtx) = s / numUpElms;
233-
});
214+
Kokkos::parallel_for(
215+
mesh.nverts(), KOKKOS_LAMBDA(const int vtx) {
216+
MeshField::Real s = 0;
217+
for (auto idx = v2e_offsets[vtx]; idx < v2e_offsets[vtx + 1]; ++idx) {
218+
const auto elm = v2e_values[idx];
219+
s += elmSize(elm);
220+
}
221+
const auto numUpElms = v2e_offsets[vtx + 1] - v2e_offsets[vtx];
222+
sizeField(vtx) = s / numUpElms;
223+
});
234224
return sizeField;
235225
}
236226

237-
template<typename EstimationT, typename OmegahMeshField, typename FieldElement>
238-
Kokkos::View<MeshField::Real*>
239-
getSprSizeField(EstimationT& e, OmegahMeshField& omf, FieldElement& coordFe) {
240-
Error errorIntegrator(e,omf);
227+
template <typename EstimationT, typename OmegahMeshField, typename FieldElement>
228+
Kokkos::View<MeshField::Real *>
229+
getSprSizeField(EstimationT &e, OmegahMeshField &omf, FieldElement &coordFe) {
230+
Error errorIntegrator(e, omf);
241231
errorIntegrator.process(coordFe);
242232
std::cout << "Error: " << errorIntegrator.r << "\n";
243233
computeSizeFactor(e, omf, coordFe, errorIntegrator);

0 commit comments

Comments
 (0)