set manifold for merged boundaries

[Weixiong Zheng]

Hi,

I just started to try manifold and tria merging functionalities. My question comes from the problem that I have to merge two curved triangulations. When refining, how do I make sure the merged boundaries are still curved?

I currently set manifold before merging triangulation but it seems not to work as desired.

See the attachment pls.

Thanks for any advice in advance.

Weixiong

[Weixiong Zheng]

I think I'd better describe a little bit more. 

I want a square domain with a cylinder inside it. What I tried doing was to get a square with a cylindrical hole merged with a cylinder. 

Here are how I did it:

1. generate the cylinder using GridGenerator::general_cell and merge_triangulation to generate the cylinder inside (has to be this way as hyperball needs to be refined and therefore can not be merged later with other triangulations)

2. set manifold on boundary of the cylinder.

3. generate the square with a cylindrical hole

4. set mainfold for it

5. merge these two triangulations

6. refine global.

Apperantly, the manifold setting is not working as what I initially thought. Any ideas how to fix it?

在 2018年3月28日星期三 UTC-7下午9:00:50，Weixiong Zheng写道：

[Weixiong Zheng]

I changed the strategy that I don't assign manifold first. I merge two triangulations. Then, I loop over cells to see if cell has a side on the merged interior boundary. If so, I set manifold id to it.

  

const Point<2> center (0,0);

const SphericalManifold<2> boundary_description(center);
tria.set_manifold (0, boundary_description);
for (typename Triangulation<2>::active_cell_iterator
       cell=tria.begin_active(); cell!=tria.end(); ++cell) {
  for (int fn=0; fn<4; ++fn) {
    auto ctr = cell->face(fn)->center();
    auto x=ctr[0],y=ctr[1];
    if (0.92<std::sqrt(x*x+y*y)<1.) {
        cell->set_all_manifold_ids(0);
        break;
    }
  }
}
tria.refine_global(2);

The result is super strange that the boundary is also exterior boundary is also curved.

在 2018年3月28日星期三 UTC-7下午9:00:50，Weixiong Zheng写道：

Hi,

[Jean-Paul Pelteret]

Hi Weixiong,

The issue that you’re seeing here is actually related to a discussion that we’re currently having regarding the behaviour of boundaries and manifolds (and will hopefully be fixed soon in the developer version). You’re hitting the issue discussed here

https://github.com/dealii/dealii/issues/5931

and we’re working through a fix for it here

https://github.com/dealii/dealii/pull/5932

The simple solution is to use a manifold ID that does not coincide with a boundary ID. Since your triangulation is not colourised, all of your boundaries will have an indicator of 0. So try to use a manifold ID of, say, 1. I think that this should fix the issue for you.

I hope this helps.

Best,

Jean-Paul

-- 
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.
<Screenshot 2018-03-28 22.08.32.png>

[Jean-Paul Pelteret]

Hi Weixiong,

So, have you looked at the FAQ as the error message suggests? There’s also actually a discussion on this point in the very first tutorial. Here’s the relevant snippet for you:

At this point, all objects created in this function will be destroyed in reverse order. Unfortunately, we defined the manifold object after the triangulation, which still has a pointer to it and the library will produce an error if the manifold object is destroyed before the triangulation. We therefore have to release it, which can be done as follows. Note that this sets the manifold object used for part "0" of the domain back to a default object, over which the triangulation has full control.

triangulation.set_manifold (0);

An alternative to doing so, and one that is frequently more convenient, would have been to declare the manifold object before the triangulation object. In that case, the triangulation would have let lose of the manifold object upon its destruction, and everything would have been fine.

Best,

Jean-Paul

On 29 Mar 2018, at 08:14, Weixiong Zheng <zwxn...@gmail.com> wrote:

Hi

It does give me straight boundary after changing to:

  

const Point<2> center (0,0);
  const SphericalManifold<2> boundary_description(center);

  //const StraightBoundary<2> boundary_description;
  //tria.set_all_manifold_ids_on_boundary(0);
  tria.set_manifold (2, boundary_description);

  
  for (typename Triangulation<2>::active_cell_iterator
       cell=tria.begin_active(); cell!=tria.end(); ++cell) {
    for (int fn=0; fn<4; ++fn) {
      auto ctr = cell->face(fn)->center();
      auto x=ctr[0],y=ctr[1];

      if (0.92<std::sqrt(x*x+y*y) && std::sqrt(x*x+y*y)<1. && !cell->at_boundary()){
        cell->set_all_manifold_ids(2);
        break;
      }
      //if (cell->at_boundary())
        //cell->set_all_manifold_ids(2);
    }
  }

But after the program finished, I got the following message for errors:

An error occurred in line <128> of file </Applications/deal.II-8.5-brew.app/Contents/Resources/brew/src/deal.II-8.5-brew/source/base/subscriptor.cc> in function

    void dealii::Subscriptor::check_no_subscribers() const

The violated condition was: 

    counter == 0

Additional information: 

    Object of class N6dealii17SphericalManifoldILi2ELi2EEE is still used by 1 other objects.

(Additional information: <none>)

See the entry in the Frequently Asked Questions of deal.II (linked to from http://www.dealii.org/) for a lot more information on what this error means and how to fix programs in which it happens.

Stacktrace:

-----------

#0  2   libdeal_II.g.8.5.0.dylib            0x000000011e04707c _ZNK6dealii11Subscriptor20check_no_subscribersEv + 3660: 2   libdeal_II.g.8.5.0.dylib            0x000000011e04707c _ZNK6dealii11Subscriptor20check_no_subscribersEv 

#1  3   libdeal_II.g.8.5.0.dylib            0x000000011e047292 _ZN6dealii11SubscriptorD2Ev + 34: 3   libdeal_II.g.8.5.0.dylib            0x000000011e047292 _ZN6dealii11SubscriptorD2Ev 

#2  4   libdeal_II.g.8.5.0.dylib            0x000000011e9c96c5 _ZN6dealii8ManifoldILi2ELi2EED2Ev + 21: 4   libdeal_II.g.8.5.0.dylib            0x000000011e9c96c5 _ZN6dealii8ManifoldILi2ELi2EED2Ev 

#3  5   step-49                             0x000000010d18f0c5 _ZN6dealii17SphericalManifoldILi2ELi2EED2Ev + 21: 5   step-49                             0x000000010d18f0c5 _ZN6dealii17SphericalManifoldILi2ELi2EED2Ev 

#4  6   step-49                             0x000000010d18eaf5 _ZN6dealii17SphericalManifoldILi2ELi2EED1Ev + 21: 6   step-49                             0x000000010d18eaf5 _ZN6dealii17SphericalManifoldILi2ELi2EED1Ev 

#5  7   step-49                             0x000000010d18d663 _Z6grid_2v + 1891: 7   step-49                             0x000000010d18d663 _Z6grid_2v 

#6  8   step-49                             0x000000010d18eb09 main + 9: 8   step-49                             0x000000010d18eb09 main 

#7  9   libdyld.dylib                       0x00007fffe1c02255 start + 1: 9   libdyld.dylib                       0x00007fffe1c02255 start 

--------------------------------------------------------

make[3]: *** [CMakeFiles/run] Abort trap: 6

make[2]: *** [CMakeFiles/run.dir/all] Error 2

make[1]: *** [CMakeFiles/run.dir/rule] Error 2

make: *** [run] Error 2

在 2018年3月28日星期三 UTC-7下午10:42:30，Jean-Paul Pelteret写道：

[Weixiong Zheng]

Both idea work fine now. Thank you so much!

在 2018年3月29日星期四 UTC-7上午11:31:23，Jean-Paul Pelteret写道：

[Jean-Paul Pelteret]

Hi Weixiong,

Great, I’m glad that you managed to fix your problem!

Best,

Jean-Paul