All,
I am having trouble with a 3D geometry/boundary description that I
have modeled after both the half hyper shell and quarter hyper shell
geometries. My geometry is a chunk of the 3D hyper shell (essentially
the 3D hyper shell cut off at z = += sqrt(3)x, z = +- sqrt(3)y) split
into 4 cells for the coarsest mesh, modeled after
GridGenerator::quarter_hyper_shell, including a colorize function. I
also implemented a boundary description based off
HalfHyperShellBoundary since I have two curved faces in my geometry
that I will want to refine and integrate over. I apply the curved
boundary description to only the two curved faces and apply no normal
flux boundary conditions on all six faces.
Unfortunately, when my ConstraintMatrix tries to close, I trigger the
exception "Cycle in constraints detected!" (line 304 of
source/lac/constraint_matrix). After doing some digging, I think the
corners are getting being overconstrained, but I can't figure out
why.
I don't have this problem in 2D, in 3D with Dirichlet boundary
conditions on the flat faces, or even in 3D with using the
quarter_hyper_shell and my own QuarterHyperShellBoundary written in
the same manner as my chunk boundary description. The first two I can
understand working, but the third case should be essentially the same
as my case, which leads me to believe the problem is in my grid or
colorize function. I can view the grid in paraview (switch to
Dirichlet on the flat faces in code below), and visually it looks as
it should, and I don't think I have vertex/cell numbering issues.
I've included a minimal code that will produce this error; a little
bit of the boundary functionality was left in for the visualization,
though it is unnecessary for producing the error. Any assistance in
finding where I'm going wrong would be greatly appreciated. Thanks in
advance,
Jennifer
-----------------------------------
#include <grid/tria.h>
#include <grid/grid_generator.h>
#include <grid/tria_boundary_lib.h>
#include <dofs/dof_handler.h>
#include <fe/fe_system.h>
#include <fe/mapping_q1.h>
#include <numerics/vectors.h>
#include <numerics/data_out.h>
#include <base/exceptions.h>
#include <base/function.h>
#include <fstream>
#include <ostream>
using namespace dealii;
template <int dim>
static void
colorize_sixty_deg_hyper_shell(Triangulation<dim> & tria,
const Point<dim>& center,
const double inner_radius,
const double outer_radius);
template<int dim>
static void sixty_deg_hyper_shell (Triangulation<dim> &tria,
const Point<dim> ¢er,
const double inner_radius,
const double outer_radius,
const unsigned int n_cells = 0,
const bool colorize = false);
template <int dim>
class SixtyDegHyperShellBoundary : public HyperShellBoundary<dim>
{
public:
SixtyDegHyperShellBoundary (const Point<dim> ¢er,
const double inner_radius,
const double outer_radius);
private:
const double inner_radius;
const double outer_radius;
};
// Implementation for 3D only
template <>
void
colorize_sixty_deg_hyper_shell(Triangulation<3> & tria,
const Point<3>& center,
const double
inner_radius,
const double
outer_radius)
{
if (tria.n_cells() != 4)
AssertThrow (false, ExcNotImplemented());
double middle = (outer_radius-inner_radius)/2e0 + inner_radius;
double eps = 1e-3*middle;
Triangulation<3>::cell_iterator cell = tria.begin();
for (;cell!=tria.end();++cell)
for(unsigned int f=0; f<GeometryInfo<3>::faces_per_cell; ++f)
{
if(!cell->face(f)->at_boundary())
continue;
double radius = cell->face(f)->center().norm() -
center.norm();
if (std::fabs(cell->face(f)->center()(2) -
sqrt(3.)*cell->face(f)->center()(0)) < eps ) // z = sqrt(3)x set
boundary 2
{
cell->face(f)->set_boundary_indicator(2);
for (unsigned int
j=0;j<GeometryInfo<3>::lines_per_face;++j)
if(cell->face(f)->line(j)->at_boundary())
if
(std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
cell->face(f)->line(j)->vertex(1).norm()) > eps)
cell->face(f)->line(j)->set_boundary_indicator(2);
}
else if (std::fabs(cell->face(f)->center()(2) +
sqrt(3.)*cell->face(f)->center()(0)) < eps) // z = -sqrt(3)x set
boundary 3
{
cell->face(f)->set_boundary_indicator(3);
for (unsigned int
j=0;j<GeometryInfo<3>::lines_per_face;++j)
if(cell->face(f)->line(j)->at_boundary())
if
(std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
cell->face(f)->line(j)->vertex(1).norm()) > eps)
cell->face(f)->line(j)->set_boundary_indicator(3);
}
else if (std::fabs(cell->face(f)->center()(2) -
sqrt(3.)*cell->face(f)->center()(1)) < eps ) // z = sqrt(3)y set
boundary 4
{
cell->face(f)->set_boundary_indicator(4);
for (unsigned int
j=0;j<GeometryInfo<3>::lines_per_face;++j)
if(cell->face(f)->line(j)->at_boundary())
if
(std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
cell->face(f)->line(j)->vertex(1).norm()) > eps)
cell->face(f)->line(j)->set_boundary_indicator(4);
}
else if (std::fabs(cell->face(f)->center()(2) +
sqrt(3.)*cell->face(f)->center()(1)) < eps ) // z = -sqrt(3)y set
boundary 5
{
cell->face(f)->set_boundary_indicator(5);
for (unsigned int
j=0;j<GeometryInfo<3>::lines_per_face;++j)
if(cell->face(f)->line(j)->at_boundary())
if
(std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
cell->face(f)->line(j)->vertex(1).norm()) > eps)
cell->face(f)->line(j)->set_boundary_indicator(5);
}
else if (radius < middle) // inner radius set boundary 0
{
cell->face(f)->set_boundary_indicator(0);
for (unsigned int
j=0;j<GeometryInfo<3>::lines_per_face;++j)
if(cell->face(f)->line(j)->at_boundary())
if
(std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
cell->face(f)->line(j)->vertex(1).norm()) < eps)
cell->face(f)->line(j)->set_boundary_indicator(0);
}
else if (radius > middle) // outer radius set boundary 1
{
cell->face(f)->set_boundary_indicator(1);
for (unsigned int
j=0;j<GeometryInfo<3>::lines_per_face;++j)
if(cell->face(f)->line(j)->at_boundary())
if
(std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
cell->face(f)->line(j)->vertex(1).norm()) < eps)
cell->face(f)->line(j)->set_boundary_indicator(1);
}
else
AssertThrow (false, ExcInternalError());
}
}
// Implementation for 3D only
template <>
void sixty_deg_hyper_shell (Triangulation<3> & tria,
const Point<3>& center,
const double inner_radius,
const double outer_radius,
const unsigned int n,
const bool colorize)
{
if (n == 0 || n == 2)
{
const double r0 = inner_radius;
const double r1 = outer_radius;
std::vector<Point<3> > vertices;
vertices.push_back (center+Point<3>( 0., 0., r0));
//0
vertices.push_back (center+Point<3>( 0., 0., r1));
//1
vertices.push_back (center+Point<3>( 0., 0.5*r0,
sqrt(3.)/2.0*r0)); //2
vertices.push_back (center+Point<3>( 0., 0.5*r1,
sqrt(3.)/2.0*r1)); //3
vertices.push_back (center+Point<3>( 0., -0.5*r0,
sqrt(3.)/2.0*r0)); //4
vertices.push_back (center+Point<3>( 0., -0.5*r1,
sqrt(3.)/2.0*r1)); //5
vertices.push_back (center+Point<3>( 0.5*r0, 0,
sqrt(3.)/2.0*r0)); //6
vertices.push_back (center+Point<3>( 0.5*r1, 0,
sqrt(3.)/2.0*r1)); //7
vertices.push_back (center+Point<3>( 1.0/sqrt(5.)*r0,
1.0/sqrt(5.)*r0, sqrt(3./5.)*r0)); //8
vertices.push_back (center+Point<3>( 1.0/sqrt(5.)*r1,
1.0/sqrt(5.)*r1, sqrt(3./5.)*r1)); //9
vertices.push_back (center+Point<3>( 1.0/sqrt(5.)*r0,
-1.0/sqrt(5.)*r0, sqrt(3./5.)*r0)); //10
vertices.push_back (center+Point<3>( 1.0/sqrt(5.)*r1,
-1.0/sqrt(5.)*r1, sqrt(3./5.)*r1)); //11
vertices.push_back (center+Point<3>( -0.5*r0, 0,
sqrt(3.)/2.0*r0)); //12
vertices.push_back (center+Point<3>( -0.5*r1, 0,
sqrt(3.)/2.0*r1)); //13
vertices.push_back (center+Point<3>( -1.0/sqrt(5.)*r0,
1.0/sqrt(5.)*r0, sqrt(3./5.)*r0)); //14
vertices.push_back (center+Point<3>( -1.0/sqrt(5.)*r1,
1.0/sqrt(5.)*r1, sqrt(3./5.)*r1)); //15
vertices.push_back (center+Point<3>( -1.0/sqrt(5.)*r0,
-1.0/sqrt(5.)*r0, sqrt(3./5.)*r0)); //16
vertices.push_back (center+Point<3>( -1.0/sqrt(5.)*r1,
-1.0/sqrt(5.)*r1, sqrt(3./5.)*r1)); //17
const int cell_vertices[4][8] = {
{12, 0, 14, 2, 13, 1, 15, 3},
{16, 4, 12, 0, 17, 5, 13, 1},
{0, 6, 2, 8, 1, 7, 3, 9},
{4, 10, 0, 6, 5, 11, 1, 7},
};
std::vector<CellData<3> > cells(4);
for (unsigned int i=0; i<4; ++i)
{
for (unsigned int j=0; j<8; ++j)
cells[i].vertices[j] = cell_vertices[i][j];
cells[i].material_id = 0;
}
tria.create_triangulation ( vertices, cells, SubCellData());
// no boundary information
}
else
{
AssertThrow(false, ExcNotImplemented());
}
if (colorize)
colorize_sixty_deg_hyper_shell(tria, center, inner_radius,
outer_radius);
}
template<int dim>
SixtyDegHyperShellBoundary<dim>::SixtyDegHyperShellBoundary (const
Point<dim> ¢er,
const double
inner_radius,
const double
outer_radius)
:
HyperShellBoundary<dim> (center),
inner_radius (inner_radius),
outer_radius (outer_radius)
{
if (dim > 2)
Assert ((inner_radius >= 0) &&
(outer_radius > 0) &&
(outer_radius > inner_radius),
ExcMessage ("Inner and outer radii must be specified
explicitly in 3d."));
}
template <int dim>
void run()
{
Triangulation<dim> triangulation;
FESystem<dim> fe(FE_Q<dim>(1), dim);
DoFHandler<dim> dof_handler (triangulation);
ConstraintMatrix constraints;
Vector<double> solution;
sixty_deg_hyper_shell (triangulation,
Point<dim>(),
0.5,
1.0,
2,
true);
static SixtyDegHyperShellBoundary<dim> boundary(Point<dim>(),
0.5,
1.0);
triangulation.set_boundary (0, boundary);
triangulation.set_boundary (1, boundary);
triangulation.refine_global(2);
dof_handler.distribute_dofs (fe);
solution.reinit(dof_handler.n_dofs());
solution = 0;
std::set<unsigned char> no_normal_flux_boundaries;
no_normal_flux_boundaries.insert (0);
VectorTools::compute_no_normal_flux_constraints
(dof_handler, 0,
no_normal_flux_boundaries,
constraints);
no_normal_flux_boundaries.insert (1);
VectorTools::compute_no_normal_flux_constraints
(dof_handler, 0,
no_normal_flux_boundaries,
constraints);
no_normal_flux_boundaries.insert (2);
VectorTools::compute_no_normal_flux_constraints
(dof_handler, 0,
no_normal_flux_boundaries,
constraints);
no_normal_flux_boundaries.insert (3);
VectorTools::compute_no_normal_flux_constraints
(dof_handler, 0,
no_normal_flux_boundaries,
constraints);
no_normal_flux_boundaries.insert (4);
VectorTools::compute_no_normal_flux_constraints
(dof_handler, 0,
no_normal_flux_boundaries,
constraints);
no_normal_flux_boundaries.insert (5);
VectorTools::compute_no_normal_flux_constraints
(dof_handler, 0,
no_normal_flux_boundaries,
constraints);
/*
VectorTools::interpolate_boundary_values (dof_handler,
2,
ZeroFunction<dim>(dim),
constraints);
VectorTools::interpolate_boundary_values (dof_handler,
3,
ZeroFunction<dim>(dim),
constraints);
VectorTools::interpolate_boundary_values (dof_handler,
4,
ZeroFunction<dim>(dim),
constraints);
VectorTools::interpolate_boundary_values (dof_handler,
5,
ZeroFunction<dim>(dim),
constraints);
*/
constraints.close();
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "solution");
data_out.build_patches();
const std::string filename = "solution.vtk";
std::ofstream output (filename.c_str());
data_out.write_vtk(output);
}
int main (int argc, char *argv[])
{
deallog.depth_console (0);
run<3>();
}
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
