domRectangle.pro

Include "domRectangle.dat";

// =================================================================================================
// FUNCTIONS
// =================================================================================================

Function{
  If (FLAG_PBM == PBM_SCATT_RECT)
    f_inc[] = Complex[Cos[WAVENUMBER*X[]], Sin[WAVENUMBER*X[]]];
    f_ref[] = AcousticFieldSoftCylinder[XYZ[]]{WAVENUMBER, R_SCA};
    f_dir[] = -f_inc[];
  EndIf
  If (FLAG_PBM == PBM_MARMOUSI)
    f_inc[] = 0;
    f_ref[] = 0;
    velocityField[] = InterpolationBilinear[ $1, $2 ]{ ListFromFile["domRectangleMarmousi.dat"] };
    cData[] = velocityField[ X[], Y[] ];
    f_dir[] = -f_inc[];
    f_sou = 1;
    For iDom In {0:N_DOM-1}
      f_sou~{iDom} = 1;
    EndFor
  EndIf
}

// =================================================================================================
// LISTS ABC/TBC
// =================================================================================================

// Is there an ABC for main domain? (1=yes, 0=no)
isEdgeABC~{0} = 1;
isEdgeABC~{1} = 1;
isEdgeABC~{2} = 1;
isEdgeABC~{3} = 1;
isEdgeRad~{0} = 0;
isEdgeRad~{1} = 0;
isEdgeRad~{2} = 0;
isEdgeRad~{3} = 0;

ListOfSubdom = {};
For xdom In {0:Nx_DOM-1}
For ydom In {0:Ny_DOM-1}
  iDom = ydom*Nx_DOM + xdom;

If (iDom % MPI_Size == MPI_Rank)
  ListOfSubdom += iDom;
EndIf

  iDomNeigh~{iDom}~{0} = (ydom-1)*Nx_DOM + xdom;  // Down
  iDomNeigh~{iDom}~{1} = ydom*Nx_DOM + (xdom+1);  // Right
  iDomNeigh~{iDom}~{2} = (ydom+1)*Nx_DOM + xdom;  // Top
  iDomNeigh~{iDom}~{3} = ydom*Nx_DOM + (xdom-1);  // Left
  iEdgeNeigh~{iDom}~{0} = 2;
  iEdgeNeigh~{iDom}~{1} = 3;
  iEdgeNeigh~{iDom}~{2} = 0;
  iEdgeNeigh~{iDom}~{3} = 1;

  // Is there an ABC or TBC? (1=yes, 0=no)
  For iEdge In {0:3}
    isEdgeABC~{iDom}~{iEdge} = 0;
    isEdgeTBC~{iDom}~{iEdge} = 1;
    isEdgeRad~{iDom}~{iEdge} = 0;
  EndFor
  If(ydom == 0)
    isEdgeABC~{iDom}~{0} = 1;
    isEdgeTBC~{iDom}~{0} = 0;
  EndIf
  If(xdom == Nx_DOM-1)
    isEdgeABC~{iDom}~{1} = 1;
    isEdgeTBC~{iDom}~{1} = 0;
  EndIf
  If(ydom == Ny_DOM-1)
    isEdgeABC~{iDom}~{2} = 1;
    isEdgeTBC~{iDom}~{2} = 0;
  EndIf
  If(xdom == 0)
    isEdgeABC~{iDom}~{3} = 1;
    isEdgeTBC~{iDom}~{3} = 0;
  EndIf

  // List of Edges for Absorbing/Transmission Boundary Condition
  ListOfEdgesWithABC~{iDom} = {};
  ListOfEdgesWithTBC~{iDom} = {};
  For iEdge In {0:3}
    If(isEdgeABC~{iDom}~{iEdge} == 1)
      ListOfEdgesWithABC~{iDom} += iEdge;
    EndIf
    If(isEdgeTBC~{iDom}~{iEdge} == 1)
      ListOfEdgesWithTBC~{iDom} += iEdge;
    EndIf
  EndFor

EndFor
EndFor

// =================================================================================================
// GROUPS
// =================================================================================================

Group {

// === REFERENCE CASE ===

  // DOMAIN
  Omega = Region[{}];

  // SOURCE
  GammaPoint = Region[{}];

  // BOUNDARY (Dirichlet BC)
  GammaD0 = Region[{}];
  GammaN0 = Region[{}];
  GammaD = Region[{}];
  GammaN = Region[{}];

  // BOUNDARY (Absorbing BC)
  For iEdge In {0:3}
    SigmaMain~{iEdge} = Region[{}];
    SigmaMainSide~{iEdge}~{0} = Region[{}];
    SigmaMainSide~{iEdge}~{1} = Region[{}];
  EndFor
  SigmaMain = Region[{}];

  // CORNER (Absorbing BC)
  CornerMain~{0} = Region[{TAG_CRN~{          0*Nx_DOM +          0 }~{0}}];
  CornerMain~{1} = Region[{TAG_CRN~{          0*Nx_DOM + (Nx_DOM-1) }~{1}}];
  CornerMain~{2} = Region[{TAG_CRN~{ (Ny_DOM-1)*Nx_DOM + (Nx_DOM-1) }~{2}}];
  CornerMain~{3} = Region[{TAG_CRN~{ (Ny_DOM-1)*Nx_DOM +          0 }~{3}}];
  CornerMain = Region[{ CornerMain~{0}, CornerMain~{1}, CornerMain~{2}, CornerMain~{3}}];

  If(FLAG_PBM == PBM_SCATT_RECT)
    GammaD += Region[{TAG_SCA}];
  EndIf

  If(FLAG_PBM == PBM_MARMOUSI)
    GammaPoint += Region[{TAG_SOU}];
  EndIf

For xdom In {0:Nx_DOM-1}
For ydom In {0:Ny_DOM-1}
  iDom = ydom*Nx_DOM + xdom;

  // DOMAIN
  Omega += Region[TAG_DOM~{iDom}];

  // BOUNDARY (Absorbing BC)
  For iEdge In {0:3}
  If(ydom == 0)
    SigmaMain~{0} += Region[{TAG_BND~{iDom}~{0}}];
    SigmaMain += Region[{TAG_BND~{iDom}~{0}}];
  EndIf
  If(xdom == Nx_DOM-1)
    SigmaMain~{1} += Region[{TAG_BND~{iDom}~{1}}];
    SigmaMain += Region[{TAG_BND~{iDom}~{1}}];
  EndIf
  If(ydom == Ny_DOM-1)
    SigmaMain~{2} += Region[{TAG_BND~{iDom}~{2}}];
    SigmaMain += Region[{TAG_BND~{iDom}~{2}}];
  EndIf
  If(xdom == 0)
    SigmaMain~{3} += Region[{TAG_BND~{iDom}~{3}}];
    SigmaMain += Region[{TAG_BND~{iDom}~{3}}];
  EndIf
  EndFor
EndFor
EndFor

  For iEdge In {0:3}
    iCornerSide1 = iEdge;
    iCornerSide2 = (iEdge+1)%4;
    SigmaMainSide~{iEdge}~{0} = Region[{CornerMain~{iCornerSide1}}];
    SigmaMainSide~{iEdge}~{1} = Region[{CornerMain~{iCornerSide2}}];
  EndFor

// === DDM CASE ===

For i0 In {0:#ListOfSubdom()-1}
iDom = ListOfSubdom(i0);

  // DOMAIN
  Omega~{iDom} = Region[TAG_DOM~{iDom}];

  // SOURCE
  GammaPoint~{iDom} = Region[{}];
  If(FLAG_PBM == PBM_MARMOUSI)
  If(iDom == Idom_SOU)
    GammaPoint~{iDom} += Region[{TAG_SOU}];
  EndIf
  EndIf

  // BOUNDARY (Dirichlet/Neumann BC)
  GammaD0~{iDom} = Region[{}];
  GammaN0~{iDom} = Region[{}];
  GammaD~{iDom} = Region[{}];
  GammaN~{iDom} = Region[{}];
  If((FLAG_PBM == PBM_SCATT_RECT) && (iDom == 0))
    GammaD~{iDom} = Region[{TAG_SCA}];
  EndIf
  TrOmegaGammaD~{iDom} = ElementsOf[ Omega~{iDom}, OnOneSideOf GammaD~{iDom} ];

  // BOUNDARY (Absorbing BC + Transmission BC)
  For iEdge In {0:3}
    iCornerSide1 = iEdge;
    iCornerSide2 = (iEdge+1)%4;
    Sigma~{iDom}~{iEdge} = Region[{TAG_BND~{iDom}~{iEdge}}];
    SigmaSide~{iDom}~{iEdge}~{0} = Region[{TAG_CRN~{iDom}~{iCornerSide1}}];
    SigmaSide~{iDom}~{iEdge}~{1} = Region[{TAG_CRN~{iDom}~{iCornerSide2}}];
    GammaD~{iDom}~{iEdge}  = Region[{}];
    GammaD0~{iDom}~{iEdge} = Region[{}];
    GammaN~{iDom}~{iEdge}  = Region[{}];
    GammaN0~{iDom}~{iEdge} = Region[{}];
  EndFor
  Sigma~{iDom} = Region[{ Sigma~{iDom}~{0}, Sigma~{iDom}~{1}, Sigma~{iDom}~{2}, Sigma~{iDom}~{3} }];

  // CORNER (Absorbing BC + Transmission BC)
  For iCorner In {0:3}
    Corner~{iDom}~{iCorner} = Region[{TAG_CRN~{iDom}~{iCorner}}];
  EndFor
  Corner~{iDom} = Region[{ Corner~{iDom}~{0}, Corner~{iDom}~{1}, Corner~{iDom}~{2}, Corner~{iDom}~{3} }];

EndFor
}