[Dune] Boundary Segments with 3D ALUGrid

[Martin Rückl]

Hi Duners,
maybe someone has had a similar problem.
Ultimately I want to achieve a rotational symmetric 3D boundary, looking 
similar to a mushroom / light bulb.
Hence, I started to create a 3D cone segment. However, my code does 
crash with an assertion error:

    spine:
    /home/mrueckl/opt/alugrid-1.52/include/serial/gitter_tetra_top.h:1116:
    ALUGridSpace::TetraTop<A>::TetraTop(int,
    ALUGridSpace::TetraTop<A>::myhface3_t*, int,
    ALUGridSpace::TetraTop<A>::myhface3_t*, int,
    ALUGridSpace::TetraTop<A>::myhface3_t*, int,
    ALUGridSpace::TetraTop<A>::myhface3_t*, int,
    ALUGridSpace::TetraTop<A>::innertetra_t*, int, double) [with A =
    ALUGridSpace::GitterBasis::Objects::TetraEmpty;
    ALUGridSpace::TetraTop<A>::myhface3_t =
    ALUGridSpace::Gitter::Geometric::hface3;
    ALUGridSpace::TetraTop<A>::innertetra_t =
    ALUGridSpace::TetraTop<ALUGridSpace::GitterBasis::Objects::TetraEmpty>]:
    Assertion `std::abs( calculatedVolume - _volume ) / _volume < 1e-10'
    failed.


In 2D I succeeded to create a segment of a circle (but only with the 
boundary segment, the boundary projection interface provided by ALUGrid 
did not create a "reasonable" grid...).
However, adding the z-component leads to the above error :-(

Any help would be appreciated!
Regards, Martin

My sample code produces the following macro grid file:

    !Tetrahedra
    4
    0.000000e+00 0.000000e+00 0.000000e+00
    1.000000e+00 0.000000e+00 0.000000e+00
    0.000000e+00 1.000000e+00 0.000000e+00
    0.000000e+00 0.000000e+00 1.000000e+00
    1
    0  1  2  3
    4
    -1  3  0  2  1
    -1  3  0  1  3
    -1  3  0  3  2
    -1  3  1  2  3
    0 -1
    1 -1
    2 -1
    3 -1

Here my sample code:

    #include <math.h>
    #include <iostream>
    #include <string>

    #include <dune/common/mpihelper.hh>
    #include <dune/common/exceptions.hh>
    #include <dune/common/fvector.hh>
    #include <dune/grid/io/file/gmshwriter.hh>

    #include <dune/grid/common/gridfactory.hh>

    #include <dune/grid/alugrid.hh>
    #include <dune/grid/alugrid/3d/alu3dgridfactory.hh>

    struct ConeBoundarySegment : public Dune::BoundarySegment<3,3>
    {

       virtual Dune::FieldVector<double,3> operator()(const
    Dune::FieldVector<double,2> &l) const
       {
         double x = l[0];    double y = l[1];

         double phi  = x*M_PI/2.0;

         Dune::FieldVector<double,3> g(0);
         g[0] = std::cos(phi) * (1.-y);
         g[1] = std::sin(phi) * (1.-y);
         g[2] = y;
         return g;
       }
    };

    int main(int argc, char** argv)
    {
       try{

       constexpr unsigned dim = 3;

       Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);

       typedef Dune::ALUGrid< dim, dim, Dune::simplex,
    Dune::nonconforming > Grid;

       Dune::GridFactory<Grid> gf;

       Dune::FieldVector<typename Grid::ctype, dim> pos;

       Dune::GeometryType type(Dune::GeometryType::simplex,dim);

        // A Tetrahedron

        // counter clockwise lower vertices
        pos[0] = 0; pos[1] =   0; pos[2] = 0; gf.insertVertex(pos);
    unsigned l0 = 0;
        pos[0] = 1; pos[1] =   0; pos[2] = 0; gf.insertVertex(pos);
    unsigned l1 = 1;
        pos[0] = 0; pos[1] =   1; pos[2] = 0; gf.insertVertex(pos);
    unsigned l2 = 2;

        // top
        pos[0] = 0;   pos[1] = 0;   pos[2] = 1; gf.insertVertex(pos); 
    unsigned t  = 3;

        gf.insertElement(type, {l0,l1,l2,t});

       // create and insert segment
       typedef Dune::shared_ptr<Dune::BoundarySegment<dim,dim> >
    SegmentPointer;
       gf.insertBoundarySegment({l1,l2,t},SegmentPointer(new
    ConeBoundarySegment()));

       Dune::shared_ptr<Grid> gridp(gf.createGrid());

       gridp->writeMacroGrid("./","macrogrid.txt");
       Dune::GmshWriter<typename Grid::LeafGridView>
    writer(gridp->leafView());

       writer.write("rr0.gmsh");
       gridp->globalRefine(1);
       writer.write("rr1.gmsh");
       gridp->globalRefine(1);
       writer.write("rr2.gmsh");
       gridp->globalRefine(1);
       writer.write("rr3.gmsh");
       gridp->globalRefine(1);
       writer.write("r4.gmsh");
       gridp->globalRefine(1);
       writer.write("r5.gmsh");
       }
       catch (Dune::Exception &e){
         std::cerr << "Dune reported error: " << e << std::endl;
       }
       catch (std::exception &e){
         std::cerr << "std::exception thrown: " << e.what()<< std::endl;
       }
       catch(...){
         std::cerr << "Unknown exception thrown!" << std::endl;
       }
    }





-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.dune-project.org/pipermail/dune/attachments/20140311/753c0a4c/attachment.htm>