4

I have a hole in an infinite plate and a code that generates a finite element mesh of it. I need more elements concentrated near the hole, because the solution varies more at this region. The issue here is that in Mesh 1 is missing some elements near the hole. The only diference between the codes is that in one i have specified MeshOrder->order, and in the other don't. Why is that?

Mesh 1 code:

Needs["NDSolve`FEM`"]
order = 2;
R = 10;
RF = 200;
mesh = ToElementMesh[
  ImplicitRegion[
   R <= Sqrt[x x + y y] && RF >= Sqrt[x x + y y] , {x, y}], {{0, 
    RF}, {0, RF}}, "MeshOrder" -> order, 
  MeshRefinementFunction -> 
   Function[{vertices, area}, 
    Block[{x, y}, {x, y} = Mean[vertices]; 
     If[(Sqrt[x x + y y] >= R && Sqrt[x x + y y] <= 1.5 R), area > 1, 
      area > 5000]]]]

mesh["Wireframe"]

enter image description here

Mesh 2 code:

mesh = ToElementMesh[
  ImplicitRegion[
   R <= Sqrt[x x + y y] && RF >= Sqrt[x x + y y] , {x, y}], {{0, 
    RF}, {0, RF}}, 
  MeshRefinementFunction -> 
   Function[{vertices, area}, 
    Block[{x, y}, {x, y} = Mean[vertices]; 
     If[(Sqrt[x x + y y] >= R && Sqrt[x x + y y] <= 1.5 R), area > 1, 
      area > 5000]]]]

mesh["Wireframe"]

enter image description here

enter image description here

Stratus
  • 2,942
  • 12
  • 24
  • What do you mean by 'discontinuities'? It's not clear what you are asking. – user21 Apr 28 '17 at 14:11
  • @user21 if you look closely in the figures you will see that in Mesh 1 are missing some elements near the hole. – Stratus Apr 28 '17 at 14:18
  • I don't see it - may you can edit you post to contain a zoomed picture that shows what you mean. – user21 Apr 28 '17 at 14:19
  • @user21 i have added a picture comparing both meshes closely. – Stratus Apr 28 '17 at 14:24
  • OK, I see. One last question - why do you specify "MeshOrder" explicitly? Any particular reason? – user21 Apr 28 '17 at 14:34
  • @user21 because i need to change every time from order=1 to order= 2 in my simulations, and this parameter (order) is called in many places in my code. – Stratus Apr 28 '17 at 14:36

1 Answers1

1

You can avoid the issue you are seeing by generating a first order boundary mesh for the numerical region and then call ToElementMesh on that. The key is to generate a first order boundary mesh - that's what ToElementMesh does in the default usage. Using "MeshOrder"->2 also generates a second order boundary mesh, which seems to interfere with the refinement function.

Needs["NDSolve`FEM`"]
order = 2;
R = 10;
RF = 200;
nr = ToNumericalRegion[
   ImplicitRegion[
    R <= Sqrt[x x + y y] && RF >= Sqrt[x x + y y], {x, y}], {{0, 
     RF}, {0, RF}}];
(* make a first order boundary mesh *)

bmesh = ToBoundaryMesh[nr, "MeshOrder" -> 1];
mesh = ToElementMesh[nr, "MeshOrder" -> order, 
   MeshRefinementFunction -> 
    Function[{vertices, area}, Block[{x, y}, {x, y} = Mean[vertices];
      If[(Sqrt[x x + y y] >= R && Sqrt[x x + y y] <= 1.5 R), area > 1,
        area > 5000]]]];
Show[mesh["Wireframe"], PlotRange -> {{-1, 20}, {-1, 20}}]

enter image description here

user21
  • 39,710
  • 8
  • 110
  • 167