Skip to content

Commit c4ac194

Browse files
roystgnrroystgnr
authored
Merge pull request #4406 from roystgnr/out_of_mesh_value_coverage
More out-of-mesh-value coverage
2 parents 21c9e74 + 1d1043a commit c4ac194

2 files changed

Lines changed: 85 additions & 32 deletions

File tree

‎src/mesh/mesh_function.C‎

Lines changed: 25 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -482,10 +482,7 @@ void MeshFunction::gradient (const Point & p,
482482

483483
const Elem * element = this->find_element(p,subdomain_ids);
484484

485-
if (!element)
486-
output.resize(0);
487-
else
488-
this->_gradient_on_elem(p, element, output);
485+
this->_gradient_on_elem(p, element, output);
489486
}
490487

491488

@@ -531,12 +528,22 @@ void MeshFunction::_gradient_on_elem (const Point & p,
531528
const Elem * element,
532529
std::vector<Gradient> & output)
533530
{
534-
libmesh_assert(element);
535-
536531
// resize the output vector to the number of output values
537532
// that the user told us
538533
output.resize (this->_system_vars.size());
539534

535+
if (!element)
536+
{
537+
libmesh_error_msg_if(!_out_of_mesh_mode,
538+
"MeshFunction couldn't find element at p=" <<
539+
p << ", but is not in out-of-mesh mode");
540+
libmesh_assert_equal_to (_out_of_mesh_value.size(),
541+
this->_system_vars.size());
542+
for (auto i : index_range(this->_system_vars))
543+
output[i] = Gradient(_out_of_mesh_value[i]);
544+
return;
545+
}
546+
540547
const unsigned int dim = element->dim();
541548

542549
// Get local coordinates to feed these into compute_data().
@@ -632,11 +639,22 @@ void MeshFunction::hessian (const Point & p,
632639
{
633640
libmesh_assert (this->initialized());
634641

642+
// resize the output vector to the number of output values
643+
// that the user told us
644+
output.resize (this->_system_vars.size());
645+
635646
const Elem * element = this->find_element(p,subdomain_ids);
636647

637648
if (!element)
638649
{
639-
output.resize(0);
650+
libmesh_error_msg_if(!_out_of_mesh_mode,
651+
"MeshFunction couldn't find element at p=" <<
652+
p << ", but is not in out-of-mesh mode");
653+
libmesh_assert_equal_to (_out_of_mesh_value.size(),
654+
this->_system_vars.size());
655+
for (auto i : index_range(this->_system_vars))
656+
output[i] = Tensor(_out_of_mesh_value[i]);
657+
return;
640658
}
641659
else
642660
{
@@ -646,11 +664,6 @@ void MeshFunction::hessian (const Point & p,
646664
<< " are not yet implemented!"
647665
<< std::endl);
648666

649-
// resize the output vector to the number of output values
650-
// that the user told us
651-
output.resize (this->_system_vars.size());
652-
653-
654667
{
655668
const unsigned int dim = element->dim();
656669

‎tests/mesh/mesh_function.C‎

Lines changed: 60 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -191,52 +191,72 @@ public:
191191

192192
EquationSystems es(mesh);
193193
System & sys = es.add_system<System>("SimpleSystem");
194-
unsigned int u_var = sys.add_variable("u", FIRST, LAGRANGE);
194+
const unsigned int u_var = sys.add_variable("u", FIRST, LAGRANGE);
195+
const unsigned int v_var = sys.add_variable("v", CONSTANT, MONOMIAL);
195196

196197
es.init();
197198
sys.project_solution(bilinear_function, nullptr, es.parameters);
198199

199-
std::vector<unsigned int> variables {u_var, libMesh::invalid_uint};
200+
std::vector<unsigned int> variables {u_var, v_var, libMesh::invalid_uint};
200201
MeshFunction mesh_function(es,
201202
*sys.current_local_solution,
202203
sys.get_dof_map(),
203204
variables);
204205
mesh_function.init();
205206

206-
DenseVector<Number> out_of_mesh_value(2);
207+
DenseVector<Number> out_of_mesh_value(3);
207208
out_of_mesh_value(0) = -99.0;
208209
out_of_mesh_value(1) = 5.25;
210+
out_of_mesh_value(1) = 42;
209211
mesh_function.enable_out_of_mesh_mode(out_of_mesh_value);
210212

211-
std::vector<Gradient> gradients;
212-
const Point p(0.4, 0.6, 0.0);
213+
// Pick an element centroid so u and v agree
214+
const Point good_p(0.375, 0.625, 0.0);
213215

214216
// On a distributed mesh not every processor may have every
215217
// element
216-
const Elem * elem = (*mesh.sub_point_locator())(p);
218+
const Elem * elem = (*mesh.sub_point_locator())(good_p);
217219
bool someone_found_elem = elem;
218220
mesh.comm().max(someone_found_elem);
219221
CPPUNIT_ASSERT(someone_found_elem);
220222

221-
mesh_function.gradient(p, 0.0, gradients);
223+
std::vector<Gradient> gradients;
224+
mesh_function.gradient(good_p, 0.0, gradients);
225+
226+
CPPUNIT_ASSERT_EQUAL(std::size_t(3), gradients.size());
222227

223228
// Let's only test our evaluation where we know we can evaluate, in parallel
224229
if (!elem || elem->processor_id() != mesh.processor_id())
225230
return;
226231

227-
CPPUNIT_ASSERT_EQUAL(std::size_t(2), gradients.size());
228-
229232
LIBMESH_ASSERT_NUMBERS_EQUAL(8.0, gradients[0](0), TOLERANCE * TOLERANCE);
230233
LIBMESH_ASSERT_NUMBERS_EQUAL(80.0, gradients[0](1), TOLERANCE * TOLERANCE);
231234
if (LIBMESH_DIM > 2)
232235
LIBMESH_ASSERT_NUMBERS_EQUAL(0.0, gradients[0](2), TOLERANCE * TOLERANCE);
233236

234-
LIBMESH_ASSERT_NUMBERS_EQUAL(out_of_mesh_value(1),
235-
gradients[1](0),
237+
// Piecewise-constant functions have zero gradients
238+
for (unsigned int d = 0; d < LIBMESH_DIM; ++d)
239+
LIBMESH_ASSERT_NUMBERS_EQUAL(0, gradients[1](d), TOLERANCE * TOLERANCE);
240+
241+
LIBMESH_ASSERT_NUMBERS_EQUAL(out_of_mesh_value(2),
242+
gradients[2](0),
236243
TOLERANCE * TOLERANCE);
237244
for (unsigned int d = 1; d < LIBMESH_DIM; ++d)
238-
LIBMESH_ASSERT_NUMBERS_EQUAL(0, gradients[1](d),
245+
LIBMESH_ASSERT_NUMBERS_EQUAL(0, gradients[2](d),
239246
TOLERANCE * TOLERANCE);
247+
248+
const Point bad_p(1.375, 0.625, 0.0);
249+
mesh_function.gradient(bad_p, 0.0, gradients);
250+
251+
for (unsigned int vn = 0; vn != 3; ++vn)
252+
{
253+
LIBMESH_ASSERT_NUMBERS_EQUAL(out_of_mesh_value(vn),
254+
gradients[vn](0),
255+
TOLERANCE * TOLERANCE);
256+
for (unsigned int d = 1; d < LIBMESH_DIM; ++d)
257+
LIBMESH_ASSERT_NUMBERS_EQUAL(0, gradients[vn](d),
258+
TOLERANCE * TOLERANCE);
259+
}
240260
}
241261

242262
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
@@ -258,20 +278,22 @@ public:
258278
EquationSystems es(mesh);
259279
System & sys = es.add_system<System>("SimpleSystem");
260280
unsigned int u_var = sys.add_variable("u", SECOND, LAGRANGE);
281+
unsigned int v_var = sys.add_variable("v", CONSTANT, MONOMIAL);
261282

262283
es.init();
263284
sys.project_solution(biquadratic_function, nullptr, es.parameters);
264285

265-
std::vector<unsigned int> variables {u_var, libMesh::invalid_uint};
286+
std::vector<unsigned int> variables {u_var, v_var, libMesh::invalid_uint};
266287
MeshFunction mesh_function(es,
267288
*sys.current_local_solution,
268289
sys.get_dof_map(),
269290
variables);
270291
mesh_function.init();
271292

272-
DenseVector<Number> out_of_mesh_value(2);
293+
DenseVector<Number> out_of_mesh_value(3);
273294
out_of_mesh_value(0) = -99.0;
274295
out_of_mesh_value(1) = 5.25;
296+
out_of_mesh_value(2) = 42;
275297
mesh_function.enable_out_of_mesh_mode(out_of_mesh_value);
276298

277299
std::vector<Tensor> hessians;
@@ -286,12 +308,12 @@ public:
286308

287309
mesh_function.hessian(p, 0.0, hessians);
288310

311+
CPPUNIT_ASSERT_EQUAL(std::size_t(3), hessians.size());
312+
289313
// Let's only test our evaluation where we know we can evaluate, in parallel
290314
if (!elem || elem->processor_id() != mesh.processor_id())
291315
return;
292316

293-
CPPUNIT_ASSERT_EQUAL(std::size_t(2), hessians.size());
294-
295317
LIBMESH_ASSERT_NUMBERS_EQUAL(15.0, hessians[0](0,0), tol);
296318
LIBMESH_ASSERT_NUMBERS_EQUAL(0.75, hessians[0](0,1), tol);
297319
LIBMESH_ASSERT_NUMBERS_EQUAL(0.75, hessians[0](1,0), tol);
@@ -301,17 +323,35 @@ public:
301323
for (unsigned int d = 0; d < LIBMESH_DIM; ++d)
302324
for (unsigned int d2 = 0; d2 < LIBMESH_DIM; ++d2)
303325
if (d>1 || d2>1)
304-
LIBMESH_ASSERT_NUMBERS_EQUAL(0, hessians[1](d,d2),
326+
LIBMESH_ASSERT_NUMBERS_EQUAL(0, hessians[0](d,d2),
305327
tol);
306328

307-
LIBMESH_ASSERT_NUMBERS_EQUAL(out_of_mesh_value(1),
308-
hessians[1](0,0),
329+
// Piecewise-constant functions have zero Hessians
330+
for (unsigned int d = 0; d < LIBMESH_DIM; ++d)
331+
for (unsigned int d2 = 0; d2 < LIBMESH_DIM; ++d2)
332+
LIBMESH_ASSERT_NUMBERS_EQUAL(0, hessians[1](d,d2), tol);
333+
334+
LIBMESH_ASSERT_NUMBERS_EQUAL(out_of_mesh_value(2),
335+
hessians[2](0,0),
309336
tol);
310337
for (unsigned int d = 0; d < LIBMESH_DIM; ++d)
311338
for (unsigned int d2 = 0; d2 < LIBMESH_DIM; ++d2)
312339
if (d || d2)
313-
LIBMESH_ASSERT_NUMBERS_EQUAL(0, hessians[1](d,d2),
340+
LIBMESH_ASSERT_NUMBERS_EQUAL(0, hessians[2](d,d2),
314341
tol);
342+
343+
const Point bad_p(1.375, 0.625, 0.0);
344+
mesh_function.hessian(bad_p, 0.0, hessians);
345+
346+
for (unsigned int vn = 0; vn != 2; ++vn)
347+
{
348+
LIBMESH_ASSERT_NUMBERS_EQUAL(out_of_mesh_value(vn),
349+
hessians[vn](0,0), tol);
350+
for (unsigned int d = 0; d < LIBMESH_DIM; ++d)
351+
for (unsigned int d2 = 0; d2 < LIBMESH_DIM; ++d2)
352+
if (d || d2)
353+
LIBMESH_ASSERT_NUMBERS_EQUAL(0, hessians[vn](d,d2), tol);
354+
}
315355
}
316356
#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
317357

0 commit comments

Comments
 (0)