OPM / IFEM

IFEM - Isogeometric Toolbox for the solution of PDEs
Other
42 stars 18 forks source link

Taylor-Hood LR basis established with wrong mesh-line multiplicity #324

Closed abdullahabdulhaque closed 5 years ago

abdullahabdulhaque commented 5 years ago

I am doing uniform refinement with LR B-splines. The polynomial degree is 2 and 3 respectively for the pressure and velocity, and their continuity is 1 everywhere except on the line from the origin to the leftmost corner, where it is 0. In the refinement, we do just 3 steps with \beta=100% (uniform). param_patch0_basis1_mesh_001 The initial mesh, which is correct. param_patch0_basis1_mesh_002 The next mesh, where the C0-line is wrong. IFEM decreases the continuity although it should not.

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<simulation>
  <geometry dim="2">
    <patchfile>Rec_Lshape.g2</patchfile>
    <!-- <raiseorder patch="1" u="1" v="1"/> -->
    <refine patch="1" type="uniform" u="7" v="7"/>
    <topologysets>
      <set name="sides" type="edge">
        <item patch="1">1 2 3 4</item>
      </set>
      <set name="horizontal" type="edge">
        <item patch="1">3 4</item>
      </set>
    </topologysets>
  </geometry>

  <stokes formulation="mixed">
    <boundaryconditions>
      <dirichlet set="sides" basis="1" comp="12" type="anasol"/>
    </boundaryconditions>
    <constrain_integrated_pressure>-0.0000052473917314548</constrain_integrated_pressure>
    <no_div_in_recovery/>
    <fluidproperties mu="1.0" rho="1.0"/>
    <anasol type="expression">
      <variables>
        rect2pol(x,y,&amp;r,&amp;t);
        r = if(below(r,1e-12),1e-12,r);
        o = 1.5*PI;
        a = 856399/1572864;
        a1 = a+1;
        a2 = a-1;
        CA = cos(a*o);
        P0 = sin(a1*t)*CA/a1-cos(a1*t)-sin(a2*t)*CA/a2+cos(a2*t);
        P1 = cos(a1*t)*CA+a1*sin(a1*t)-cos(a2*t)*CA-a2*sin(a2*t);
        P2 = -a1*sin(a1*t)*CA+a1*a1*cos(a1*t)+a2*sin(a2*t)*CA-a2*a2*cos(a2*t);
        P3 = -a1*a1*(cos(a1*t)*CA+a1*sin(a1*t))+a2*a2*(cos(a2*t)*CA+a2*sin(a2*t));
        u = (a1*sin(t)*P0+cos(t)*P1)*pow(r,a);
        ux = (a*cos(2*t)*P1+(a1*a2*P0-P2)*0.5*sin(2*t))*pow(r,a2);
        uy = (pow(cos(t),2)*P2+a*sin(2*t)*P1-a1*P0*(a2*pow(cos(t),2)-a))*pow(r,a2);
        v = (sin(t)*P1-a1*cos(t)*P0)*pow(r,a);
        vx = (-pow(sin(t),2)*P2+a*sin(2*t)*P1-a1*P0*(a2*pow(cos(t),2)+1))*pow(r,a2);
        vy = -ux;
        p = pow(r,a-1)*(a1*a1*P1+P3)/a2;
      </variables>
      <primary>u|v</primary>
      <scalarprimary>p</scalarprimary>
      <secondary>ux|uy|vx|vy</secondary>
    </anasol>
    <source type="expression">
            rect2pol(x,y,&amp;r,&amp;t);
            r = if(below(r,1e-16),1e-16,r);
            o = 1.5*PI;
            a = 856399/1572864;
            a1 = a+1;
            a2 = a-1;
            CA = cos(a*o);
            P0 = sin(a1*t)*CA/a1-cos(a1*t)-sin(a2*t)*CA/a2+cos(a2*t);
            P1 = cos(a1*t)*CA+a1*sin(a1*t)-cos(a2*t)*CA-a2*sin(a2*t);
            P2 = -a1*sin(a1*t)*CA+a1*a1*cos(a1*t)+a2*sin(a2*t)*CA-a2*a2*cos(a2*t);
            P3 = -a1*a1*(cos(a1*t)*CA+a1*sin(a1*t))+a2*a2*(cos(a2*t)*CA+a2*sin(a2*t));
            P4 = pow(a1,3)*(sin(a1*t)*CA-a1*cos(a1*t))+pow(a2,3)*(a2*cos(a2*t)-sin(a2*t)*CA);
            S1 = a1*a1*P1+P3;
            S2 = a1*a1*P2+P4;
            D2u = (cos(t)*P3+a2*sin(t)*P2+a1*a1*(cos(t)*P1+a2*sin(t)*P0))*pow(r,a-2);
            D2v = (sin(t)*P3-a2*cos(t)*P2+a1*a1*(sin(t)*P1-a2*cos(t)*P0))*pow(r,a-2);
            px = (S1*cos(t)-S2*sin(t)/a2)*pow(r,a-2);
            py = (S1*sin(t)+S2*cos(t)/a2)*pow(r,a-2);
      -D2u+px | -D2v+py
    </source>
  </stokes>

  <discretization>
    <nGauss>5</nGauss>
  </discretization>

  <postprocessing>
    <projection>
      <CGL2/>
    </projection>
  </postprocessing>

  <linearsolver>
    <rtol>1e-10</rtol>
    <ilu_fill_level>1</ilu_fill_level>
  </linearsolver>

  <adaptive>
    <maxstep value="3"/>
    <maxdof value="40000"/>
    <beta type="symmetrized" value="100"/>
    <errtol value="0.00000001"/>
    <use_norm value="1"/>
    <scheme value="isotropic_function"/>
    <scheme>isotropic_function</scheme>
    <store_eps_mesh/>
    <store_errors/>
  </adaptive>

</simulation>

200 1 0 0
2 0
5 3
0 0 0 1 1 2 2 2
3 3
0 0 0 1 1 1
0 -1
0 -0.5
0 0
0.5 0
1 0
-0.5 -1
-0.5 -0.25
-0.5 0.5
0.25 0.5
1 0.5
-1 -1
-1 0
-1 1
0 1
1 1
``
akva2 commented 5 years ago

now fixed.