[CIG-SEISMO] Negative Jacobian in specfem2d

[Dimitri Komatitsch]

Hi again Jami,

PS: the two codes (Gmesh and SPECFEM2D) have a slightly different way of 
defining the Jacobian in the case of high order and highly distorted 
(curved, or very distorted) elements. Thus you can get close to zero in 
both, but slightly above zero in one and slightly below in the other; 
but in both cases the mesh is unusable because it contains a few 
elements that are too flat and will make the simulation blow up anyway.

Please let us know,
Thanks,
Best wishes,

Dimitri.

On 03/08/2016 15:43, Dimitri Komatitsch wrote:
>
> Hi Jami,
>
> Thanks! Probably mesh quality that is not good enough, two elements that
> are almost flat. Try to see the quality histogram in Gmsh to see if
> these two elements are poor (in which case even if you manage to run
> that mesh in SPECFEM2D it will likely blow up). There is also a tool for
> that in SPECFEM2D:
>
> src/auxiliaries/check_quality_external_mesh.f90
>
> If it turns out to be a bug in the code (unlikely, but worth checking)
> please let us know.
>
> Thanks,
> Best wishes,
>
> Dimitri.
>
> On 03/08/2016 10:01, Jami Johnson wrote:
>> Hello,
>> I have created a very dense mesh with Gmsh, and when I check for the
>> values of the Jacobian in Gmsh they are all positive.  However, when I
>> run it in specfem2d, an error is thrown due to a negative Jacobian.  It
>> seems to be a very small percentage of my elements causing a problem.  I
>> believe all of my line loops are oriented correctly.  Are there any
>> other potential causes of a negative Jacobian?
>>
>> Error in specfem2d:
>> --------------------------
>>  2  elements have a negative Jacobian, out of      1313871
>>
>> geo file:
>> -----------
>> lc=0.00003;
>> xmax=0.015;
>> xmin=0.0;
>> ymin=-0.015;
>> ymax=0.0;
>> Point(1) = {xmax,ymax , 0, lc};
>> Point(2) = {xmax, ymin, 0, lc};
>> Point(3) = {xmin, ymin, 0, lc};
>> Point(4) = {xmin, ymax, 0, lc};
>> Line(1) = {1, 4};
>> Line(2) = {4, 3};
>> Line(3) = {3, 2};
>> Line(4) = {2, 1};
>> R=0.00025;
>> Centx=0.0075;
>> Centy=-0.005;
>> Point(5) = {Centx, R+Centy, 0, lc};
>> Point(6) = {R+Centx, Centy, 0, lc};
>> Point(7) = {-R+Centx, Centy, 0, lc};
>> Point(8) = {Centx, -R+Centy, 0, lc};
>> Point(9) = {Centx, Centy, 0, lc};
>> Circle(5) = {5, 9, 7};
>> Circle(6) = {7, 9, 8};
>> Circle(7) = {8, 9, 6};
>> Circle(8) = {6, 9, 5};
>> Line Loop(9) = {5, 6, 7, 8};
>> Line Loop(10) = {1, 2, 3, 4};
>>
>> Plane Surface(10) = {10, 9};
>> Plane Surface(11) = {9};
>>
>> Recombine Surface{10,11};
>> Mesh.SubdivisionAlgorithm=1;
>> Physical Line("Top") = {1};
>> Physical Line("Left") = {2};
>> Physical Line("Bottom") = {3};
>> Physical Line("Right") = {4};
>> Physical Surface("M1") = {10};
>> Physical Surface("M2") = {11};
>>
>>
>>
>>
>> _______________________________________________
>> CIG-SEISMO mailing list
>> CIG-SEISMO at geodynamics.org
>> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-seismo
>>
>

-- 
Dimitri Komatitsch
CNRS Research Director (DR CNRS), Laboratory of Mechanics and Acoustics,
UPR 7051, Marseille, France    http://komatitsch.free.fr