Skip to content

Commit 2c17963

Browse files
author
Andrea Mola
committed
hypersingular and stringly singular quadrature added for gradients computations. detached from toolkit
1 parent 82ad3d2 commit 2c17963

14 files changed

Lines changed: 720 additions & 1016 deletions

‎CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -124,6 +124,8 @@ if(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/Doxyfile.in)
124124
endif()
125125

126126

127+
message(STATUS "DEAL_II_INCLUDE_DIRECTORIES=${deal.II_INCLUDE_DIRECTORIES}")
128+
127129
ADD_CUSTOM_TARGET(indent
128130
WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}
129131
COMMAND ./scripts/indent

‎include/ass_leg_function.h‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,6 @@
2020
#include <deal.II/lac/precondition.h>
2121
#include <deal.II/lac/solver_control.h>
2222
#include <deal.II/lac/solver_gmres.h>
23-
#include <deal.II/lac/sparse_matrix.h>
2423
#include <deal.II/lac/vector.h>
2524
// #include <deal.II/grid/tria_boundary_lib.h>
2625

‎include/bem_problem.h‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -346,10 +346,17 @@ class BEMProblem : public ParameterAcceptor
346346

347347
TrilinosWrappers::MPI::Vector system_rhs;
348348

349+
/// solution and alpha vectors
349350
TrilinosWrappers::MPI::Vector sol;
350351
TrilinosWrappers::MPI::Vector alpha;
352+
/// an alternatively computed alpha vector (obtained with geometric computations)
351353
TrilinosWrappers::MPI::Vector hyp_alpha;
352-
std::vector<TrilinosWrappers::MPI::Vector> C_ii;
354+
/// a set of distributed vectors which contain all the entries of the
355+
/// C_ij tensor appearing in the hypersingular BIE
356+
std::vector<TrilinosWrappers::MPI::Vector> C_ij;
357+
/// a set of distributed vectors which contain all the entries of the
358+
/// b_i vector appearing in the hypersingular BIE
359+
std::vector<TrilinosWrappers::MPI::Vector> b_i;
353360

354361
mutable TrilinosWrappers::MPI::Vector serv_phi;
355362
mutable TrilinosWrappers::MPI::Vector serv_dphi_dn;

‎include/singular_kernel_integral.h‎

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -62,7 +62,7 @@ SingularKernelIntegral(typename DoFHandler<dim-1, dim>::active_cell_iterator in_
6262
Mapping<dim-1, dim> &in_mapping,
6363
Point<dim-1> &in_eta);
6464

65-
Tensor<1,dim> evaluate_Vk_integrals();
65+
6666

6767
Tensor<1,dim> evaluate_free_term_b();
6868

@@ -73,8 +73,6 @@ std::vector<Tensor<1,dim> > evaluate_VkNj_integrals();
7373
std::vector<Tensor<1,dim> > evaluate_WkNj_integrals();
7474

7575

76-
double evaluate_integral();
77-
7876
private:
7977
const typename DoFHandler<dim-1, dim>::active_cell_iterator &cell;
8078
FiniteElement<dim-1, dim> &fe;

‎source/bem_problem.cc‎

Lines changed: 174 additions & 136 deletions
Large diffs are not rendered by default.

‎source/computational_domain.cc‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -521,7 +521,7 @@ ComputationalDomain<dim>::create_initial_mesh()
521521

522522
GridTools::delete_unused_vertices(vertices, cells, subcelldata);
523523
GridTools::consistently_order_cells(cells);
524-
tria.create_triangulation_compatibility(vertices, cells, subcelldata);
524+
tria.create_triangulation(vertices, cells, subcelldata);
525525

526526
static const Point<dim> center = Point<dim>();
527527
static const SphericalManifold<dim - 1, dim> manifold(center);

‎source/singular_kernel_integral.cc‎

Lines changed: 38 additions & 874 deletions
Large diffs are not rendered by default.

‎tests/hypersingular_integral.cc‎

Lines changed: 121 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,121 @@
1+
//-----------------------------------------------------------
2+
//
3+
// Copyright (C) 2014 by the deal.II authors
4+
//
5+
// This file is subject to LGPL and may not be distributed
6+
// without copyright and license information. Please refer
7+
// to the file deal.II/doc/license.html for the text and
8+
// further information on this license.
9+
//
10+
//-----------------------------------------------------------
11+
12+
13+
#include "singular_kernel_integral.h"
14+
#include "tests.h"
15+
16+
int
17+
main()
18+
{
19+
initlog();
20+
21+
//we start selecting a mapping degree
22+
double mapping_degree = 1;
23+
24+
// and now we generate a grid and use it to distribute dofs
25+
Triangulation<2,3> triangulation;
26+
27+
std::vector<Point<3> > vertices;
28+
std::vector<CellData<2> > cells;
29+
SubCellData subcelldata;
30+
31+
// Flat
32+
33+
vertices.resize(4);
34+
vertices[0](0)=-1.0;
35+
vertices[0](1)=-1.0;
36+
vertices[0](2)=0.0;
37+
vertices[1](0)=1.5;
38+
vertices[1](1)=-1.0;
39+
vertices[1](2)=0.0;
40+
vertices[2](0)=-1.0;
41+
vertices[2](1)=1.0;
42+
vertices[2](2)=0.0;
43+
vertices[3](0)=0.5;
44+
vertices[3](1)=1.0;
45+
vertices[3](2)=0.0;
46+
47+
48+
cells.resize(1);
49+
50+
cells[0].vertices[0]=0;
51+
cells[0].vertices[1]=1;
52+
cells[0].vertices[2]=2;
53+
cells[0].vertices[3]=3;
54+
55+
GridTools::delete_unused_vertices (vertices, cells, subcelldata);
56+
GridTools::consistently_order_cells(cells);
57+
58+
// here's the triangulation set up
59+
triangulation.create_triangulation(vertices, cells, subcelldata );
60+
61+
62+
// we will need a codimension one finite element and
63+
// the corresponding dof handler
64+
// we use zeroth order fe_dgq because in the integral the hypersingular kernel must not be multiplied by the shape function. Thus,
65+
// we multiply it by a constant shape function
66+
FE_DGQ<2,3> fe(0);
67+
DoFHandler<2,3> dof_handler(triangulation);
68+
dof_handler.distribute_dofs(fe);
69+
MappingQ<2,3> mapping(mapping_degree);
70+
71+
// The test cases are taken from
72+
//M. Guiggiani, G. Krishnasamy, T. J. Rudolphi & F. J. Rizzo,A general algorithm for the numerical solution of hypersingular boundary integral equations, JOURNAL OF APPLIED MECHANICS-TRANSACTIONS OF THE ASME,vol. 59,pp 604-614,tot.pag 11,1992
73+
74+
// Case 4.1a
75+
Point<2> singularity_location_in_parametric_plane_1(0.5,0.5);
76+
//4.1b
77+
Point<2> singularity_location_in_parametric_plane_2(0.66*.5+.5,0.5);
78+
//4.1c
79+
Point<2> singularity_location_in_parametric_plane_3(0.885764071856287,0.66*.5+.5);
80+
81+
std::vector<Tensor<1,3> > I;
82+
83+
for (const auto &cell : dof_handler.active_cell_iterators())
84+
{
85+
deallog<<"Case: 4.1a"<<std::endl;
86+
SingularKernelIntegral sing_kernel_integrator_1(cell, fe, mapping, singularity_location_in_parametric_plane_1);
87+
I = sing_kernel_integrator_1.evaluate_VkNj_integrals();
88+
//std::vector<Tensor<1,3> > I = sing_kernel_integrator.evaluate_WkNj_integrals();
89+
deallog<<std::setprecision(15);
90+
deallog<<"Result : "<<I[0][2]*4*dealii::numbers::PI<<" Exact: "<<-5.749236751228080<<std::endl;
91+
deallog<<"Abs Err: "<<-5.749236751228080-I[0][2]*4*dealii::numbers::PI<<std::endl;
92+
deallog<<"Rel Err: "<<fabs(-5.749236751228080-I[0][2]*4*dealii::numbers::PI)/fabs(-5.749236751228080)<<std::endl;
93+
94+
95+
deallog<<"Case: 4.1b"<<std::endl;
96+
SingularKernelIntegral sing_kernel_integrator_2(cell, fe, mapping, singularity_location_in_parametric_plane_2);
97+
I = sing_kernel_integrator_2.evaluate_VkNj_integrals();
98+
//std::vector<Tensor<1,3> > I = sing_kernel_integrator.evaluate_WkNj_integrals();
99+
deallog<<std::setprecision(15);
100+
deallog<<"Result : "<<I[0][2]*4*dealii::numbers::PI<<" Exact: "<<-9.154585469918885<<std::endl;
101+
deallog<<"Abs Err: "<<-9.154585469918885-I[0][2]*4*dealii::numbers::PI<<std::endl;
102+
deallog<<"Rel Err: "<<fabs(-9.154585469918885-I[0][2]*4*dealii::numbers::PI)/fabs(-9.154585469918885)<<std::endl;
103+
104+
105+
deallog<<"Case: 4.1c"<<std::endl;
106+
SingularKernelIntegral sing_kernel_integrator_3(cell, fe, mapping, singularity_location_in_parametric_plane_3);
107+
I = sing_kernel_integrator_3.evaluate_VkNj_integrals();
108+
//std::vector<Tensor<1,3> > I = sing_kernel_integrator.evaluate_WkNj_integrals();
109+
deallog<<std::setprecision(15);
110+
deallog<<"Result : "<<I[0][2]*4*dealii::numbers::PI<<" Exact: "<<-15.32849545090306<<std::endl;
111+
deallog<<"Abs Err: "<<-15.32849545090306-I[0][2]*4*dealii::numbers::PI<<std::endl;
112+
deallog<<"Rel Err: "<<fabs(-15.32849545090306-I[0][2]*4*dealii::numbers::PI)/fabs(-15.32849545090306)<<std::endl;
113+
114+
115+
116+
117+
}
118+
119+
120+
121+
}

‎tests/hypersingular_integral.output‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,13 @@
1+
2+
DEAL::Case: 4.1a
3+
DEAL::Result : -5.74923675122754 Exact: -5.74923675122808
4+
DEAL::Abs Err: -5.42677014436777e-13
5+
DEAL::Rel Err: 9.43911405841578e-14
6+
DEAL::Case: 4.1b
7+
DEAL::Result : -9.15458546967812 Exact: -9.15458546991889
8+
DEAL::Abs Err: -2.40769182369149e-10
9+
DEAL::Rel Err: 2.63003915535329e-11
10+
DEAL::Case: 4.1c
11+
DEAL::Result : -15.3284954556581 Exact: -15.3284954509031
12+
DEAL::Abs Err: 4.75505679276012e-09
13+
DEAL::Rel Err: 3.10210275234806e-10

‎tests/hypersingular_integral_curved_cell.cc‎

Lines changed: 154 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,154 @@
1+
//-----------------------------------------------------------
2+
//
3+
// Copyright (C) 2014 by the deal.II authors
4+
//
5+
// This file is subject to LGPL and may not be distributed
6+
// without copyright and license information. Please refer
7+
// to the file deal.II/doc/license.html for the text and
8+
// further information on this license.
9+
//
10+
//-----------------------------------------------------------
11+
12+
#include <deal.II/opencascade/manifold_lib.h>
13+
#include <deal.II/opencascade/utilities.h>
14+
#include "singular_kernel_integral.h"
15+
#include "tests.h"
16+
17+
int
18+
main()
19+
{
20+
initlog();
21+
22+
//we start selecting a mapping degree
23+
double mapping_degree = 2;
24+
//we also import a CAD geometry of a portion of a cylinder to integrate on curved cell
25+
std::string cad_filename = std::string(SOURCE_DIR)+"/testData/Revolution_1.iges";
26+
TopoDS_Shape surface = OpenCASCADE::read_IGES(cad_filename, 1e-3);
27+
OpenCASCADE::NormalToMeshProjectionManifold<2,3> manifold(surface, 1e-7);
28+
// the next lines are in case we want to run the test with NURBSPatchManifold
29+
// and --- possibly --- with MappingManifold. But there is some work to do to make it work
30+
// std::vector<TopoDS_Face> occ_faces;
31+
// std::vector<TopoDS_Edge> occ_edges;
32+
// std::vector<TopoDS_Vertex> occ_vertices;
33+
// OpenCASCADE::extract_geometrical_shapes(surface, occ_faces, occ_edges, occ_vertices);
34+
// OpenCASCADE::NURBSPatchManifold<2,3> manifold(occ_faces[0], 1e-7);
35+
36+
37+
// and now we generate a grid and use it to distribute dofs
38+
Triangulation<2,3> triangulation;
39+
40+
std::vector<Point<3> > vertices;
41+
std::vector<CellData<2> > cells;
42+
SubCellData subcelldata;
43+
44+
// Flat
45+
46+
vertices.resize(4);
47+
vertices[0](0)=1.0;
48+
vertices[0](1)=0.0;
49+
vertices[0](2)=0.0;
50+
vertices[1](0)=1.0;
51+
vertices[1](1)=2.0;
52+
vertices[1](2)=0.0;
53+
vertices[2](0)=0.0;
54+
vertices[2](1)=0.0;
55+
vertices[2](2)=1.0;
56+
vertices[3](0)=0.0;
57+
vertices[3](1)=2.0;
58+
vertices[3](2)=1.0;
59+
60+
61+
cells.resize(1);
62+
63+
cells[0].vertices[0]=0;
64+
cells[0].vertices[1]=1;
65+
cells[0].vertices[2]=2;
66+
cells[0].vertices[3]=3;
67+
68+
GridTools::delete_unused_vertices (vertices, cells, subcelldata);
69+
GridTools::consistently_order_cells(cells);
70+
71+
// here's the triangulation set up
72+
triangulation.create_triangulation(vertices, cells, subcelldata );
73+
74+
triangulation.set_all_manifold_ids(0);
75+
triangulation.set_manifold(0, manifold);
76+
77+
78+
79+
// we will need a codimension one finite element and
80+
// the corresponding dof handler
81+
// we use zeroth order fe_dgq because in the integral the hypersingular kernel must not be multiplied by the shape function. Thus,
82+
// we multiply it by a constant shape function
83+
FE_DGQ<2,3> fe(0);
84+
DoFHandler<2,3> dof_handler(triangulation);
85+
dof_handler.distribute_dofs(fe);
86+
MappingQ<2,3> mapping(mapping_degree);
87+
88+
// The test cases are taken from
89+
//M. Guiggiani, G. Krishnasamy, T. J. Rudolphi & F. J. Rizzo,A general algorithm for the numerical solution of hypersingular boundary integral equations, JOURNAL OF APPLIED MECHANICS-TRANSACTIONS OF THE ASME,vol. 59,pp 604-614,tot.pag 11,1992
90+
91+
// Case 4.2a
92+
Point<2> singularity_location_in_parametric_plane_1(0.5,0.5);
93+
// Case 4.2b
94+
Point<2> singularity_location_in_parametric_plane_2(0.66*.5+.5,0.0*.5+.5);
95+
// Case 4.2c
96+
Point<2> singularity_location_in_parametric_plane_3(0.66*.5+.5,0.66*.5+.5);
97+
98+
99+
std::vector<Tensor<1,3> > I;
100+
101+
for (const auto &cell : dof_handler.active_cell_iterators())
102+
{
103+
deallog<<"Case: 4.1a"<<std::endl;
104+
SingularKernelIntegral sing_kernel_integrator_1(cell, fe, mapping, singularity_location_in_parametric_plane_1);
105+
I = sing_kernel_integrator_1.evaluate_VkNj_integrals();
106+
//std::vector<Tensor<1,3> > I = sing_kernel_integrator.evaluate_WkNj_integrals();
107+
deallog<<std::setprecision(15);
108+
deallog<<"Result : "<<I[0][2]*4*dealii::numbers::PI<<" Exact: "<<-0.343807*4*dealii::numbers::PI<<std::endl;
109+
deallog<<"Abs Err: "<<-0.343807*4*dealii::numbers::PI-I[0][2]*4*dealii::numbers::PI<<std::endl;
110+
deallog<<"Rel Err: "<<fabs(-0.343807*4*dealii::numbers::PI-I[0][2]*4*dealii::numbers::PI)/fabs(-0.343807*4*dealii::numbers::PI)<<std::endl;
111+
112+
113+
deallog<<"Case: 4.1b"<<std::endl;
114+
SingularKernelIntegral sing_kernel_integrator_2(cell, fe, mapping, singularity_location_in_parametric_plane_2);
115+
I = sing_kernel_integrator_2.evaluate_VkNj_integrals();
116+
//std::vector<Tensor<1,3> > I = sing_kernel_integrator.evaluate_WkNj_integrals();
117+
deallog<<std::setprecision(15);
118+
deallog<<"Result : "<<I[0][2]*4*dealii::numbers::PI<<" Exact: "<<-0.497099*4*dealii::numbers::PI<<std::endl;
119+
deallog<<"Abs Err: "<<-0.497099*4*dealii::numbers::PI-I[0][2]*4*dealii::numbers::PI<<std::endl;
120+
deallog<<"Rel Err: "<<fabs(-0.497099*4*dealii::numbers::PI-I[0][2]*4*dealii::numbers::PI)/fabs(-0.497099*4*dealii::numbers::PI)<<std::endl;
121+
122+
123+
deallog<<"Case: 4.1c"<<std::endl;
124+
SingularKernelIntegral sing_kernel_integrator_3(cell, fe, mapping, singularity_location_in_parametric_plane_3);
125+
I = sing_kernel_integrator_3.evaluate_VkNj_integrals();
126+
//std::vector<Tensor<1,3> > I = sing_kernel_integrator.evaluate_WkNj_integrals();
127+
deallog<<std::setprecision(15);
128+
deallog<<"Result : "<<I[0][2]*4*dealii::numbers::PI<<" Exact: "<<-0.877214*4*dealii::numbers::PI<<std::endl;
129+
deallog<<"Abs Err: "<<-0.877214*4*dealii::numbers::PI-I[0][2]*4*dealii::numbers::PI<<std::endl;
130+
deallog<<"Rel Err: "<<fabs(-0.877214*4*dealii::numbers::PI-I[0][2]*4*dealii::numbers::PI)/fabs(-0.877214*4*dealii::numbers::PI)<<std::endl;
131+
132+
133+
134+
// strong integrals case a
135+
// deallog<<"Abs Err: "<<-2.114175-I[0][0]*4*dealii::numbers::PI<<std::endl;
136+
// deallog<<"Rel Err: "<<fabs(-2.114175-I[0][0]*4*dealii::numbers::PI)/fabs(-2.114175)<<std::endl;
137+
138+
139+
// strong integrals case b
140+
// deallog<<"Abs Err: "<<-1.935711-I[0][0]*4*dealii::numbers::PI<<std::endl;
141+
// deallog<<"Rel Err: "<<fabs(-1.935711-I[0][0]*4*dealii::numbers::PI)/fabs(-1.935711)<<std::endl;
142+
143+
// strong integrals case c
144+
// deallog<<"Abs Err: "<<.8790179-I[0][0]*4*dealii::numbers::PI<<std::endl;
145+
// deallog<<"Rel Err: "<<fabs(.8790179-I[0][0]*4*dealii::numbers::PI)/fabs(.8790179)<<std::endl;
146+
147+
148+
149+
150+
}
151+
152+
153+
154+
}

0 commit comments

Comments
 (0)