Merge pull request #93 from WenjieYao/emscatter

Update electromagnetic scattering tutorial

Author: fverdugo

.github/workflows/ci.yml

@@ -93,7 +93,7 @@ jobs:
             ${{ runner.os }}-test-
             ${{ runner.os }}-
       - uses: julia-actions/julia-buildpkg@v1
-      - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl darcy.jl  inc_navier_stokes.jl stokes.jl poisson_dev_fe.jl 
+      - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl darcy.jl  inc_navier_stokes.jl stokes.jl poisson_dev_fe.jl emscatter.jl
   tutorials_set_4:
     name: Tutorials4 ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
     runs-on: ${{ matrix.os }}

Project.toml

@@ -8,6 +8,7 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
+GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8"
 IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -16,9 +17,11 @@ NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
 Gridap = "0.16.3"
+GridapGmsh = "0.4"
 julia = "1.3"
 
 [extras]

assets/emscatter/Illustration.png

@@ -0,0 +1,78 @@
+using Gmsh
+import Gmsh: gmsh
+function MeshGenerator(L,H,xc,r,hs,d_pml,lc)
+    gmsh.initialize()
+    gmsh.option.setNumber("General.Terminal", 1)
+    gmsh.option.setNumber("Mesh.Algorithm", 6)
+    gmsh.clear()
+    gmsh.model.add("geometry")
+
+    # Add points
+    gmsh.model.geo.addPoint(-L/2-d_pml, -H/2-d_pml, 0, lc, 1)
+    gmsh.model.geo.addPoint( L/2+d_pml, -H/2-d_pml, 0, lc, 2)
+    gmsh.model.geo.addPoint( L/2+d_pml,  hs       , 0, lc, 3)
+    gmsh.model.geo.addPoint(-L/2-d_pml,  hs       , 0, lc, 4)
+    gmsh.model.geo.addPoint( L/2+d_pml,  H/2+d_pml, 0, lc, 5)
+    gmsh.model.geo.addPoint(-L/2-d_pml,  H/2+d_pml, 0, lc, 6)
+    gmsh.model.geo.addPoint( xc[1]-r  ,  xc[2]    , 0, lc, 7)
+    gmsh.model.geo.addPoint( xc[1]    ,  xc[2]    , 0, lc, 8)
+    gmsh.model.geo.addPoint( xc[1]+r  ,  xc[2]    , 0, lc, 9)
+    # Add lines
+    gmsh.model.geo.addLine( 1,  2,  1)
+    gmsh.model.geo.addLine( 2,  3,  2)
+    gmsh.model.geo.addLine( 3,  4,  3)
+    gmsh.model.geo.addLine( 1,  4,  4)
+    gmsh.model.geo.addLine( 3,  5,  5)
+    gmsh.model.geo.addLine( 5,  6,  6)
+    gmsh.model.geo.addLine( 4,  6,  7)
+    gmsh.model.geo.addCircleArc( 7, 8, 9, 8)
+    gmsh.model.geo.addCircleArc( 9, 8, 7, 9)
+    # Construct curve loops and surfaces 
+    gmsh.model.geo.addCurveLoop([1, 2, 3, -4], 1)
+    gmsh.model.geo.addCurveLoop([5, 6,-7, -3], 2)
+    gmsh.model.geo.addCurveLoop([8, 9], 3)
+    gmsh.model.geo.addPlaneSurface([1,3], 1)
+    gmsh.model.geo.addPlaneSurface([2], 2)
+    gmsh.model.geo.addPlaneSurface([3], 3)
+    # Physical groups
+    #gmsh.model.addPhysicalGroup(0, [1,2,3,4,5,6], 1)
+    #gmsh.model.setPhysicalName(0, 1, "DirichletNodes")
+    gmsh.model.addPhysicalGroup(1, [1,6], 2)
+    gmsh.model.setPhysicalName(1, 2, "DirichletEdges")
+    gmsh.model.addPhysicalGroup(0, [7,9], 3)
+    gmsh.model.setPhysicalName(0, 3, "CylinderNodes")
+    gmsh.model.addPhysicalGroup(1, [8,9], 4)
+    gmsh.model.setPhysicalName(1, 4, "CylinderEdges")
+    gmsh.model.addPhysicalGroup(2, [3], 5)
+    gmsh.model.setPhysicalName(2, 5, "Cylinder")
+    gmsh.model.addPhysicalGroup(2, [1,2], 6)
+    gmsh.model.setPhysicalName(2, 7, "Air")
+    gmsh.model.addPhysicalGroup(1, [3], 7)
+    gmsh.model.setPhysicalName(1, 7, "Source")
+    gmsh.model.geo.synchronize()
+    
+    gmsh.model.mesh.setPeriodic(1, [2], [4],
+                            [1, 0, 0, L+2*d_pml, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1])
+    gmsh.model.mesh.setPeriodic(1, [5], [7],
+                            [1, 0, 0, L+2*d_pml, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1])
+    
+    # We can then generate a 2D mesh...
+    gmsh.model.mesh.generate(2)
+    # ... and save it to disk
+    gmsh.write("geometry.msh")
+    gmsh.finalize()
+end
+
+# Geometry parameters
+λ = 1.0          # Wavelength (arbitrary unit)
+L = 4.0          # Width of the area
+H = 6.0          # Height of the area
+xc = [0 -1.0]    # Center of the cylinder
+r = 1.0          # Radius of the cylinder
+hs = 2.0         # y-position of the source (plane wave)
+d_pml = 0.8      # Thickness of the PML
+
+resol = 20.0      # Number of points per wavelength
+lc = λ/resol      # Characteristic length
+
+MeshGenerator(L,H,xc,r,hs,d_pml,lc)

assets/emscatter/MeshGenerator.jl

@@ -16,6 +16,7 @@ files = [
   "Stokes equation" => "stokes.jl",
   "Isotropic damage model" => "isotropic_damage.jl",
   "Fluid-Structure Interaction"=>"fsi_tutorial.jl",
+  "Electromagnetic scattering in 2D"=>"emscatter.jl",
   "Low-level API Poisson equation"=>"poisson_dev_fe.jl",
   "On using DrWatson.jl"=>"validation_DrWatson.jl"]