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
