domRectangle.geo

Include "domRectangle.dat";

Mesh.CharacteristicLengthMin = LC;
Mesh.CharacteristicLengthMax = LC;

// =================================================================================================
// Build SCATTERER or SOURCE
// =================================================================================================

If(FLAG_PBM == PBM_SCATT_RECT)
  p0 = newp; Point(p0) = {     0,     0, 0};
  p1 = newp; Point(p1) = { R_SCA,     0, 0};
  p2 = newp; Point(p2) = {     0, R_SCA, 0};
  p3 = newp; Point(p3) = {-R_SCA,     0, 0};
  p4 = newp; Point(p4) = {     0,-R_SCA, 0};
  c1 = newl; Circle(c1) = {p1, p0, p2};
  c2 = newl; Circle(c2) = {p2, p0, p3};
  c3 = newl; Circle(c3) = {p3, p0, p4};
  c4 = newl; Circle(c4) = {p4, p0, p1};
  llSca = newll; Line loop(llSca) = {-c1, -c2, -c3, -c4};
EndIf

If(FLAG_PBM == PBM_MARMOUSI)
  p0 = newp; Point(p0) = {X_SOU, Y_SOU, 0};
  Physical Point(TAG_SOU) = p0;
EndIf

// =================================================================================================
// Build POINTS, STRAIGHT LINES, RECTANGLES
// =================================================================================================

For i In {0:Nx_DOM}
For j In {0:Ny_DOM}
  p~{i}~{j} = newp;
  If((i == 1) && (j == 1))
    xAdd = TWIST;
    yAdd = 0;
  ElseIf((i == 1) && (j == 2))
    xAdd = 0;
    yAdd = TWIST;
  ElseIf((i == 2) && (j == 2))
    xAdd = -TWIST;
    yAdd = 0;
  ElseIf((i == 2) && (j == 1))
    xAdd = 0;
    yAdd = -TWIST;
  Else
    xAdd = 0;
    yAdd = 0;
  EndIf
  Point(p~{i}~{j}) = {X0_DOM + Lx*i/Nx_DOM + xAdd, Y0_DOM + Ly*j/Ny_DOM + yAdd, 0};
EndFor
EndFor

For i In {0:Nx_DOM}
For j In {0:Ny_DOM}
  If (i < Nx_DOM)
    lx~{i}~{j} = newl;
    Line(lx~{i}~{j}) = {p~{i}~{j}, p~{i+1}~{j}};
  EndIf
  If (j < Ny_DOM)
    ly~{i}~{j} = newl;
    Line(ly~{i}~{j}) = {p~{i}~{j}, p~{i}~{j+1}};
  EndIf
EndFor
EndFor

For i In {0:Nx_DOM-1}
For j In {0:Ny_DOM-1}
  ll~{i}~{j} = newll;
  Line loop(ll~{i}~{j}) = {lx~{i}~{j}, ly~{i+1}~{j}, -lx~{i}~{j+1}, -ly~{i}~{j}};
  s~{i}~{j} = news;
  If ((i == 0) && (j == 0) && (FLAG_PBM == PBM_SCATT_RECT))
    Plane Surface(s~{i}~{j}) = {ll~{i}~{j}, llSca};
  Else
    Plane Surface(s~{i}~{j}) = {ll~{i}~{j}};
  EndIf
  If (FLAG_PBM == PBM_MARMOUSI)
  If ((i == I_SOU) && (j == J_SOU))
    Point {p0} In Surface {s~{i}~{j}};
  EndIf
  EndIf
EndFor
EndFor

// =================================================================================================
// Generate MESH for each subdomain and reference domain
// =================================================================================================

If(StrCmp(OnelabAction, "check"))
  Mesh 2;
EndIf

For i In {0:Nx_DOM-1}
For j In {0:Ny_DOM-1}
  idom = j*Nx_DOM + i;
  Delete Physicals;

  Physical Point(TAG_CRN~{idom}~{0}) = p~{i}~{j};
  Physical Point(TAG_CRN~{idom}~{1}) = p~{i+1}~{j};
  Physical Point(TAG_CRN~{idom}~{2}) = p~{i+1}~{j+1};
  Physical Point(TAG_CRN~{idom}~{3}) = p~{i}~{j+1};
  Physical Line(TAG_BND~{idom}~{0}) = lx~{i}~{j};
  Physical Line(TAG_BND~{idom}~{1}) = ly~{i+1}~{j};
  Physical Line(TAG_BND~{idom}~{2}) = lx~{i}~{j+1};
  Physical Line(TAG_BND~{idom}~{3}) = ly~{i}~{j};
  Physical Surface(TAG_DOM~{idom}) = s~{i}~{j};

  If((i == 0) && (j == 0) && (FLAG_PBM == PBM_SCATT_RECT))
    Physical Line(TAG_SCA) = {-c1, -c2, -c3, -c4};
  EndIf
  If(FLAG_PBM == PBM_MARMOUSI)
  If((i == I_SOU) && (j == J_SOU))
    Physical Point(TAG_SOU) = p0;
  EndIf
  EndIf

  Printf("Saving subdomain %g...", idom);
  Save StrCat(MSH_NAME, Sprintf("_%g.msh", idom));
  Printf("Done.");
EndFor
EndFor

Delete Physicals;

If(FLAG_PBM == PBM_SCATT_RECT)
  Physical Line(TAG_SCA) = {-c1, -c2, -c3, -c4};
EndIf
If(FLAG_PBM == PBM_MARMOUSI)
  Physical Point(TAG_SOU) = p0;
EndIf

For i In {0:Nx_DOM-1}
For j In {0:Ny_DOM-1}
  idom = j*Nx_DOM + i;
  Physical Surface(TAG_DOM~{idom}) = s~{i}~{j};

  If(j == 0)
    Physical Line(TAG_BND~{idom}~{0}) = lx~{i}~{j};
  EndIf
  If(i == Nx_DOM-1)
    Physical Line(TAG_BND~{idom}~{1}) = ly~{i+1}~{j};
  EndIf
  If(j == Ny_DOM-1)
    Physical Line(TAG_BND~{idom}~{2}) = lx~{i}~{j+1};
  EndIf
  If(i == 0)
    Physical Line(TAG_BND~{idom}~{3}) = ly~{i}~{j};
  EndIf

EndFor
EndFor

Physical Point(TAG_CRN~{          0*Nx_DOM +          0 }~{0}) = p~{0     }~{0     };
Physical Point(TAG_CRN~{          0*Nx_DOM + (Nx_DOM-1) }~{1}) = p~{Nx_DOM}~{0     };
Physical Point(TAG_CRN~{ (Ny_DOM-1)*Nx_DOM + (Nx_DOM-1) }~{2}) = p~{Nx_DOM}~{Ny_DOM};
Physical Point(TAG_CRN~{ (Ny_DOM-1)*Nx_DOM +          0 }~{3}) = p~{0     }~{Ny_DOM};

Printf("Saving domain ...");
Save StrCat(MSH_NAME, Sprintf(".msh"));
Printf("Done.");