/
mapping_manifold_07.cc
96 lines (75 loc) · 2.72 KB
/
mapping_manifold_07.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
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 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.
//
// ---------------------------------------------------------------------
// Try to compute the area of a circle/sphere using JxW values.
#include <deal.II/base/utilities.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_manifold.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria.h>
#include "../tests.h"
template <int dim, int spacedim>
void
test()
{
std::ostream &out = deallog.get_file_stream();
out << "# dim=" << dim << ", spacedim=" << spacedim << std::endl;
Triangulation<dim, spacedim> triangulation;
Point<spacedim> center;
center[0] = 1.5;
center[1] = 2.5;
double radius = 1.0;
static const PolarManifold<dim, spacedim> manifold(center);
GridGenerator::hyper_ball(triangulation, center, radius);
triangulation.set_all_manifold_ids_on_boundary(0);
triangulation.set_manifold(0, manifold);
for (unsigned int cycle = 0; cycle < 4; ++cycle)
{
MappingManifold<dim, spacedim> map_manifold;
FE_Q<dim, spacedim> fe(1);
const QGauss<dim - 1> quad(3);
FEFaceValues<dim, spacedim> fe_v(map_manifold,
fe,
quad,
update_JxW_values);
double area = 0;
for (typename Triangulation<dim, spacedim>::active_cell_iterator cell =
triangulation.begin_active();
cell != triangulation.end();
++cell)
for (const unsigned int f : GeometryInfo<dim>::face_indices())
if (cell->face(f)->at_boundary())
{
fe_v.reinit(cell, f);
for (unsigned int i = 0; i < quad.size(); ++i)
area += fe_v.JxW(i);
}
deallog << "Cycle : " << cycle << std::endl;
deallog << "Surface Area : " << area << std::endl;
deallog << "Error : " << (area - (dim - 1) * 2 * numbers::PI)
<< std::endl;
triangulation.refine_global(1);
}
}
int
main()
{
initlog();
test<2, 2>();
test<3, 3>();
return 0;
}