solve on mesh imported through opencascade interface

[thomas stephens]

I am trying to combine techniques from tutorial 54 and tutorial 38 - namely, I would like to solve a surface diffusion problem on a geometry that I bring in through the opencascade interface.  

I have created a geometry (using Trelis/Cubit) and have the mesh .ucd file, as well as the .iges file (attached as sphere_mesh.ucd and sphere.iges).  I can build a GridIn<2,3> object from sphere_mesh.ucd and attach that to a Triangulation, and I can create a TopoDS_Shape object from sphere.iges, and I can attach a NormalProjectionBoundary<2,3> normal_projector to the triangulation - at least to the extent that I can then call triangulation.refine_global(n) and obtain a refinement of sphere_mesh.ucd that continues to yield better approximations to a sphere.  (This all occurs in the make_grid_and_dofs() function in laplace_on_surface.cc)

The problem I am having shows up in the assemble_system() function in laplace_on_surface.cc.  The error I get is 

 

An error occurred in line <53> of file </home/tds3/software/dealii/dealii-8.4.1/source/grid/tria_boundary.cc> in function 
void dealii::Boundary<dim, spacedim>::get_intermediate_points_on_quad(const typename dealii::Triangulation<dim, spacedim>::quad_iterator&, std::vector<dealii::Point<spacedim> >&) const [with int dim = 2; int spacedim = 3; typename dealii::Triangulation<dim, spacedim>::quad_iterator = dealii::TriaIterator<dealii::CellAccessor<2, 3> >]
The violated condition was: 
    false
The name and call sequence of the exception was:
    ExcPureFunctionCalled()
You (or a place in the library) are trying to call a function that is declared as a virtual function in a base class but that has not been overridden in your derived class.

The entire stacktrace is:

#0  /usr/local/dealii/lib/libdeal_II.g.so.8.4.1: dealii::Boundary<2, 3>::get_intermediate_points_on_quad(dealii::TriaIterator<dealii::CellAccessor<2, 3> > const&, std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > >&) const
#1  /usr/local/dealii/lib/libdeal_II.g.so.8.4.1: 
#2  /usr/local/dealii/lib/libdeal_II.g.so.8.4.1: dealii::MappingQGeneric<2, 3>::add_quad_support_points(dealii::TriaIterator<dealii::CellAccessor<2, 3> > const&, std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > >&) const
#3  /usr/local/dealii/lib/libdeal_II.g.so.8.4.1: dealii::MappingQGeneric<2, 3>::compute_mapping_support_points(dealii::TriaIterator<dealii::CellAccessor<2, 3> > const&) const
#4  /usr/local/dealii/lib/libdeal_II.g.so.8.4.1: dealii::MappingQGeneric<2, 3>::fill_fe_values(dealii::TriaIterator<dealii::CellAccessor<2, 3> > const&, dealii::CellSimilarity::Similarity, dealii::Quadrature<2> const&, dealii::Mapping<2, 3>::InternalDataBase const&, dealii::internal::FEValues::MappingRelatedData<2, 3>&) const
#5  /usr/local/dealii/lib/libdeal_II.g.so.8.4.1: dealii::MappingQ<2, 3>::fill_fe_values(dealii::TriaIterator<dealii::CellAccessor<2, 3> > const&, dealii::CellSimilarity::Similarity, dealii::Quadrature<2> const&, dealii::Mapping<2, 3>::InternalDataBase const&, dealii::internal::FEValues::MappingRelatedData<2, 3>&) const
#6  /usr/local/dealii/lib/libdeal_II.g.so.8.4.1: dealii::FEValues<2, 3>::do_reinit()
#7  /usr/local/dealii/lib/libdeal_II.g.so.8.4.1: void dealii::FEValues<2, 3>::reinit<dealii::DoFHandler, false>(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 3>, false> > const&)
#8  ./laplace_on_surface: laplace_on_surface::LaplaceBeltramiProblem<3>::assemble_system()
#9  ./laplace_on_surface: laplace_on_surface::LaplaceBeltramiProblem<3>::run()
#10  ./laplace_on_surface: main

It seems to me that I am missing something in the way I bring in the geometry back in make_grid_and_dofs() because I can generate a spherical mesh using the GridGenerator objects (as is done in tutorial step 38, and everything works fine.  My suspicion is that I am not setting manifold ids correctly in my current implementation above.  

Also, please suggest appropriate tags for this question - I thought this should get an opencascade tag, but that's not a 'suggested tag'.  

Thank you, 

Tom

[luca.heltai]

Dear Tom,

you’re not doing anything wrong. The Manifold/Boundary interfaces are undergoing a lot of restructuring, and you hit a problem of backward compatibility implementation…

Your issue is in the following lines of code, that were inserted to maintain backward compatibility between Manifold and Boundary objects:

https://github.com/dealii/dealii/pull/2418/files#diff-c743c51b4f2c95b842d938d7f90da390L3545

The above pull request is work in progress, as it breaks a lot of backward compatibilities. I hope to have some time soon to work on it...

In particular, you see that if an object is derived by a Boundary class, then a different code is executed w.r.t pure Manifold classes (Boundary is derived from Manifold, but knows of being a codimension one surface).

Unfortunately the OpenCascade wrappers were derived from Boundary (they were introduced before the Manifold concept was formulated), but they do not implement all function calls

get_intermediate_points_on_line/quad/hex

which is what triggers your exception. In your case it should be sufficient to copy the class you are interested in and make it derived from FlatManifold instead of
Boundary. Everything should then work out of the box.

By the way, you are hitting this problem because you try to use higher order mappings on co-dimension one surfaces, defined through iges files. This should works, but its highly untested… you may hit some high order mapping bugs we are trying to address… I hope things will be resolved soon.

Can you open an issue, so that I can refer back to it when things are resolved?

Best,

Luca.

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <CMakeLists.txt><laplace_on_surface.cc><sphere_mesh.ucd><sphere.iges>

[thomas stephens]

Luca,

Thank you for getting back to me so quickly.  Unfortunately I am still stuck.

I sort of understand what you are saying here: 

 it should be sufficient to copy the class you are interested in and make it derived from FlatManifold instead of Boundary

but I don't know exactly which class you are referring to or how to implement that.   My guess is something like:

  static OpenCASCADE::NormalProjectionBoundary<2,3> normal_projector(cad_surface, tolerance);

 FlatManifold<2,3> normal_projector_as_FlatManifold = FlatManifold(normal_projector_copy);

 triangulation.set_manifold(1,normal_projector_as_FlatManifold);

where I am confused as to how to obtain normal_projector_copy.  I'm also uncertain if this is the class(?) you are referring to.

Additionally, you provided the helpful remark: 

By the way, you are hitting this problem because you try to use higher order mappings on co-dimension one surfaces, defined through iges files

I see that I am using the default degree=2 to define the order of my mappings on the codimension one surface.  I ran the laplace_on_surface code using degree = 1 and this moved my error to another location in the code base!  This is a positive sign, but doesn't get me all the way home.  I've attached a modified laplace_on_surface.cc file which has the following two features:

1. an option in make_grid_and_dofs() which specifies where to get the mesh and manifold description from.  The string manifold_method specifies "cad" or "GridGenerator", and I have implemented a similar codimension one sphere from the GridGenerator utility.  
2. I have supplied the degree argument in the instantiation of the LaplaceProblem<spacedim> class in the main() function, this currently is set to 1 (rather than the default, degree=2).  

When I run the new function on the sphere generated using GridGenerator life is good (for degree 1 and 2).  When I run the new function with the sphere generated by  manifold_method="cad" and degree 2 I get the same old problem (as expected) - when I run with degree 1, the old problem appears to go away but the solver does not converge. 

----------------------------------------------------
Exception on processing: 
--------------------------------------------------------
An error occurred in line <606> of file </usr/local/dealii/include/deal.II/lac/solver_cg.h> in function
    void dealii::SolverCG<VectorType>::solve(const MatrixType&, VectorType&, const VectorType&, const PreconditionerType&) [with MatrixType = dealii::SparseMatrix<double>; PreconditionerType = dealii::PreconditionSSOR<>; VectorType = dealii::Vector<double>]

The violated condition was: 
    false
The name and call sequence of the exception was:

    SolverControl::NoConvergence (it, res)
Additional Information: 
Iterative method reported convergence failure in step 3074. The residual in the last step was 2.09928.

This is unexpected because the sphere generated from cad should be basically the same as the sphere generated by GridGenerator.  It seems to me that I still cannot treat meshes and manifolds brought in through cad as though they were the same as ones generated by GridGenerator.  

Thank you for your time,

Tom

You received this message because you are subscribed to a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/XG6J8oX1XxQ/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.