add topology optimization tutorial

.github/workflows/ci.yml

@@ -124,7 +124,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 isotropic_damage.jl validation_DrWatson.jl interpolation_fe.jl poisson_transient.jl
+      - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl isotropic_damage.jl validation_DrWatson.jl interpolation_fe.jl poisson_transient.jl TopOptEMFocus.jl
   tutorials_mpi:
     name: TutorialsMPI ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
     runs-on: ${{ matrix.os }}

Project.toml

@@ -4,12 +4,15 @@ authors = ["Santiago Badia <santiago.badia@monash.edu>", "Francesc Verdugo <fver
 version = "0.17.0"
 
 [deps]
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 GridapDistributed = "f9701e48-63b3-45aa-9a63-9bc6c271f355"
 GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8"
+GridapMakie = "41f30b06-6382-4b60-a5f7-79d86b35bf5d"
 GridapP4est = "c2c8e14b-f5fd-423d-9666-1dd9ad120af9"
 GridapPETSc = "bcdc36c2-0c3e-11ea-095a-c9dadae499f1"
 IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
@@ -18,12 +21,15 @@ LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
+NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
 Gridap = "0.17.9"

assets/TopOptEMFocus/MeshGenerator.jl

@@ -0,0 +1,106 @@
+using Gmsh
+import Gmsh: gmsh
+
+struct RecCirGeometry
+    L::Float64           # Length of the normal region  
+    h1::Float64          # Height of normal region
+    h2::Float64          # Height of the region above source
+    rt::Float64          # Radius of the target location
+    rd::Float64          # Radius of design domain
+    rs::Float64          # Radius of smallest distance circle
+    dpml::Float64        # Thickness of the PML
+    # Characteristic length (controls the resolution, smaller the finer)
+    l1::Float64          # Normal region
+    l2::Float64          # Design domain
+end
+
+
+function MeshGenerator(geo_param::RecCirGeometry, meshfile_name::String)
+    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(-geo_param.L/2-geo_param.dpml, -geo_param.dpml, 0, geo_param.l1, 1)
+    gmsh.model.geo.addPoint( geo_param.L/2+geo_param.dpml, -geo_param.dpml, 0, geo_param.l1, 2)
+    gmsh.model.geo.addPoint( geo_param.L/2+geo_param.dpml,  geo_param.h1, 0, geo_param.l1, 3)
+    gmsh.model.geo.addPoint( geo_param.L/2+geo_param.dpml,  geo_param.h1+geo_param.h2+geo_param.dpml, 0, geo_param.l1, 4)
+    gmsh.model.geo.addPoint(-geo_param.L/2-geo_param.dpml,  geo_param.h1+geo_param.h2+geo_param.dpml, 0, geo_param.l1, 5)
+    gmsh.model.geo.addPoint(-geo_param.L/2-geo_param.dpml,  geo_param.h1, 0, geo_param.l1, 6)
+    gmsh.model.geo.addPoint(0,  geo_param.h1/2, 0, geo_param.l2, 7)
+    gmsh.model.geo.addPoint(-geo_param.rs,  geo_param.h1/2, 0, geo_param.l2, 8)
+    gmsh.model.geo.addPoint( geo_param.rs,  geo_param.h1/2, 0, geo_param.l2, 9)
+    gmsh.model.geo.addPoint(-geo_param.rd,  geo_param.h1/2, 0, geo_param.l2, 10)
+    gmsh.model.geo.addPoint( geo_param.rd,  geo_param.h1/2, 0, geo_param.l2, 11)
+    # 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( 4,  5,  4)
+    gmsh.model.geo.addLine( 6,  5,  5)
+    gmsh.model.geo.addLine( 1,  6,  6)
+    gmsh.model.geo.addLine( 6,  3,  7)
+    gmsh.model.geo.addCircleArc( 8, 7, 9, 8)
+    gmsh.model.geo.addCircleArc( 9, 7, 8, 9)
+    gmsh.model.geo.addCircleArc(10, 7,11,10)
+    gmsh.model.geo.addCircleArc(11, 7,10,11)
+    # Construct curve loops and surfaces 
+    gmsh.model.geo.addCurveLoop([8, 9], 1)
+    gmsh.model.geo.addPlaneSurface([1], 1)
+    gmsh.model.geo.addCurveLoop([10, 11], 2)
+    gmsh.model.geo.addPlaneSurface([2, 1], 2)
+    gmsh.model.geo.addCurveLoop([1, 2, -7,-6], 4)
+    gmsh.model.geo.addPlaneSurface([4, 2], 4)
+    gmsh.model.geo.addCurveLoop([3, 4,-5, 7], 5)
+    gmsh.model.geo.addPlaneSurface([5], 5)
+    # Physical groups
+    gmsh.model.addPhysicalGroup(0, [1, 2, 4, 5], 1)
+    gmsh.model.setPhysicalName(0, 1, "DirichletNodes")
+    gmsh.model.addPhysicalGroup(1, [1, 4], 2)
+    gmsh.model.setPhysicalName(1, 2, "DirichletEdges")
+    gmsh.model.addPhysicalGroup(0, [8, 9, 10, 11], 3)
+    gmsh.model.setPhysicalName(0, 3, "DesignNodes")
+    gmsh.model.addPhysicalGroup(1, [8, 9, 10, 11], 4)
+    gmsh.model.setPhysicalName(1, 4, "DesignEdges")
+    gmsh.model.addPhysicalGroup(2, [2], 5)
+    gmsh.model.setPhysicalName(2, 5, "Design")
+    gmsh.model.addPhysicalGroup(2, [1], 6)
+    gmsh.model.setPhysicalName(2, 6, "Center")
+    gmsh.model.addPhysicalGroup(2, [3,4,5], 7)
+    gmsh.model.setPhysicalName(2, 7, "Air")
+    gmsh.model.addPhysicalGroup(1, [7], 9)
+    gmsh.model.setPhysicalName(1, 9, "Source")
+    gmsh.model.geo.synchronize()
+
+    # Set periodic mesh on the left and right side
+    gmsh.model.mesh.setPeriodic(1, [2], [6],
+                            [1, 0, 0, geo_param.L+2*geo_param.dpml, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1])
+    gmsh.model.mesh.setPeriodic(1, [3], [5],
+                            [1, 0, 0, geo_param.L+2*geo_param.dpml, 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(meshfile_name)
+    gmsh.finalize()
+end
+
+L = 600
+h1 = 600
+h2 = 200
+rt = 150
+rd = 100
+rs = 10
+
+dpml = 300
+
+l1 = 20
+l2 = 1
+
+meshfile = "RecCirGeometry.msh"
+geo_param = RecCirGeometry(L, h1, h2, rt, rd, rs, dpml, l1, l2)
+MeshGenerator(geo_param, meshfile)

deps/build.jl

@@ -21,7 +21,8 @@ files = [
   "On using DrWatson.jl"=>"validation_DrWatson.jl",
   "Interpolation of CellFields"=>"interpolation_fe.jl",
   "Poisson equation on parallel distributed-memory computers"=>"poisson_distributed.jl",
-  "Transient Poisson equation"=>"poisson_transient.jl"]
+  "Transient Poisson equation"=>"poisson_transient.jl",
+  "Topology optimization"=>"TopOptEMFocus.jl"]
 
 Sys.rm(notebooks_dir;recursive=true,force=true)
 for (i,(title,filename)) in enumerate(files)