-
Notifications
You must be signed in to change notification settings - Fork 737
/
find_cell_4.cc
74 lines (60 loc) · 2.08 KB
/
find_cell_4.cc
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
// ---------------------------------------------------------------------
//
// Copyright (C) 2003 - 2020 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------
// take a 3d mesh and check that we can find an arbitrary point's cell
// in it. this presently fails since the point sits right on the edge
// of the domain, but for different reasons that find_cell_5
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include "../tests.h"
void
check(Triangulation<3> &tria)
{
Point<3> p(0.75, 0, 0);
Triangulation<3>::active_cell_iterator cell =
GridTools::find_active_cell_around_point(tria, p);
deallog << cell << std::endl;
for (const unsigned int v : GeometryInfo<3>::vertex_indices())
deallog << "<" << cell->vertex(v) << "> ";
deallog << std::endl;
// Transform back and forth
Point<3> pp = StaticMappingQ1<3>::mapping.transform_unit_to_real_cell(
cell,
GeometryInfo<3>::project_to_unit_cell(
StaticMappingQ1<3>::mapping.transform_real_to_unit_cell(cell, p)));
AssertThrow(p.distance(pp) < 1e-15, ExcInternalError());
}
int
main()
{
initlog();
try
{
Triangulation<3> coarse_grid;
GridGenerator::hyper_cube(coarse_grid);
coarse_grid.refine_global(3);
check(coarse_grid);
}
catch (const std::exception &exc)
{
// we shouldn't get here...
deallog << "Caught an error..." << std::endl;
deallog << exc.what() << std::endl;
}
}