diff --git a/data/compass.geo b/data/compass.geo new file mode 100644 index 00000000..90572084 --- /dev/null +++ b/data/compass.geo @@ -0,0 +1,118 @@ +SetFactory("OpenCASCADE"); + +order = 1; + +R = 1; +r = 0.2; + +Point(1) = {0,0,0}; + +Point(2) = {r/Sqrt(2),r/Sqrt(2),0}; +Point(3) = {-r/Sqrt(2),r/Sqrt(2),0}; +Point(4) = {-r/Sqrt(2),-r/Sqrt(2),0}; +Point(5) = {r/Sqrt(2),-r/Sqrt(2),0}; + +Point(6) = {R,0,0}; +Point(7) = {R/Sqrt(2),R/Sqrt(2),0}; +Point(8) = {0,R,0}; +Point(9) = {-R/Sqrt(2),R/Sqrt(2),0}; +Point(10) = {-R,0,0}; +Point(11) = {-R/Sqrt(2),-R/Sqrt(2),0}; +Point(12) = {0,-R,0}; +Point(13) = {R/Sqrt(2),-R/Sqrt(2),0}; + +Line(1) = {1,2}; +Line(2) = {1,3}; +Line(3) = {1,4}; +Line(4) = {1,5}; + +Line(5) = {1,6}; +Line(6) = {1,8}; +Line(7) = {1,10}; +Line(8) = {1,12}; + +Line(9) = {2,6}; +Line(10) = {2,8}; +Line(11) = {3,8}; +Line(12) = {3,10}; +Line(13) = {4,10}; +Line(14) = {4,12}; +Line(15) = {5,12}; +Line(16) = {5,6}; + +Line(17) = {6,7}; +Line(18) = {7,8}; +Line(19) = {8,9}; +Line(20) = {9,10}; +Line(21) = {10,11}; +Line(22) = {11,12}; +Line(23) = {12,13}; +Line(24) = {13,6}; + +Transfinite Curve{1:24} = 2; + +Physical Curve("ENE") = {17}; +Physical Curve("NNE") = {18}; +Physical Curve("NNW") = {19}; +Physical Curve("WNW") = {20}; +Physical Curve("WSW") = {21}; +Physical Curve("SSW") = {22}; +Physical Curve("SSE") = {23}; +Physical Curve("ESE") = {24}; + +Curve Loop(1) = {9,17,18,-10}; +Curve Loop(2) = {11,19,20,-12}; +Curve Loop(3) = {13,21,22,-14}; +Curve Loop(4) = {15,23,24,-16}; + +Plane Surface(1) = {1}; +Plane Surface(2) = {2}; +Plane Surface(3) = {3}; +Plane Surface(4) = {4}; + +Transfinite Surface{1} = {2,6,7,8}; +Transfinite Surface{2} = {3,8,9,10}; +Transfinite Surface{3} = {4,10,11,12}; +Transfinite Surface{4} = {5,12,13,6}; +Recombine Surface{1:4}; + +Physical Surface("Base") = {1,2,3,4}; + +Curve Loop(5) = {1,10,-6}; +Plane Surface(5) = {5}; +Physical Surface("N Even") = {5}; + +Curve Loop(6) = {6,-11,-2}; +Plane Surface(6) = {6}; +Physical Surface("N Odd") = {6}; + +Curve Loop(7) = {2,12,-7}; +Plane Surface(7) = {7}; +Physical Surface("W Even") = {7}; + +Curve Loop(8) = {7,-13,-3}; +Plane Surface(8) = {8}; +Physical Surface("W Odd") = {8}; + +Curve Loop(9) = {3,14,-8}; +Plane Surface(9) = {9}; +Physical Surface("S Even") = {9}; + +Curve Loop(10) = {8,-15,-4}; +Plane Surface(10) = {10}; +Physical Surface("S Odd") = {10}; + +Curve Loop(11) = {4,16,-5}; +Plane Surface(11) = {11}; +Physical Surface("E Even") = {11}; + +Curve Loop(12) = {5,-9,-1}; +Plane Surface(12) = {12}; +Physical Surface("E Odd") = {12}; + +// Generate 2D mesh +Mesh 2; +SetOrder order; +Mesh.MshFileVersion = 2.2; + +Save "compass.msh"; diff --git a/data/compass.mesh b/data/compass.mesh new file mode 100644 index 00000000..8bff51ce --- /dev/null +++ b/data/compass.mesh @@ -0,0 +1,96 @@ +MFEM mesh v1.3 + +# +# MFEM Geometry Types (see mesh/geom.hpp): +# +# POINT = 0 +# SEGMENT = 1 +# TRIANGLE = 2 +# SQUARE = 3 +# TETRAHEDRON = 4 +# CUBE = 5 +# PRISM = 6 +# + +dimension +2 + +elements +12 +10 2 7 0 1 +11 2 0 7 2 +12 2 9 0 2 +13 2 0 9 3 +14 2 11 0 3 +15 2 0 11 4 +16 2 5 0 4 +17 2 0 5 1 +9 3 1 5 6 7 +9 3 2 7 8 9 +9 3 3 9 10 11 +9 3 4 11 12 5 + +attribute_sets +16 +"Base" 1 9 +"E Even" 1 16 +"E Odd" 1 17 +"East" 2 16 17 +"N Even" 1 10 +"N Odd" 1 11 +"North" 2 10 11 +"Rose" 8 10 11 12 13 14 15 16 17 +"Rose Even" 4 10 12 14 16 +"Rose Odd" 4 11 13 15 17 +"S Even" 1 14 +"S Odd" 1 15 +"South" 2 14 15 +"W Even" 1 12 +"W Odd" 1 13 +"West" 2 12 13 + +boundary +8 +1 1 5 6 +2 1 6 7 +3 1 7 8 +4 1 8 9 +5 1 9 10 +6 1 10 11 +7 1 11 12 +8 1 12 5 + +bdr_attribute_sets +13 +"Boundary" 8 1 2 3 4 5 6 7 8 +"ENE" 1 1 +"ESE" 1 8 +"Eastern Boundary" 2 1 8 +"NNE" 1 2 +"NNW" 1 3 +"Northern Boundary" 2 2 3 +"SSE" 1 7 +"SSW" 1 6 +"Southern Boundary" 2 6 7 +"WNW" 1 4 +"WSW" 1 5 +"Western Boundary" 2 4 5 + +vertices +13 +2 +0 0 +0.14142136 0.14142136 +-0.14142136 0.14142136 +-0.14142136 -0.14142136 +0.14142136 -0.14142136 +1 0 +0.70710678 0.70710678 +0 1 +-0.70710678 0.70710678 +-1 0 +-0.70710678 -0.70710678 +0 -1 +0.70710678 -0.70710678 + +mfem_mesh_end diff --git a/data/compass.msh b/data/compass.msh new file mode 100644 index 00000000..f9229490 --- /dev/null +++ b/data/compass.msh @@ -0,0 +1,62 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$PhysicalNames +17 +1 1 "ENE" +1 2 "NNE" +1 3 "NNW" +1 4 "WNW" +1 5 "WSW" +1 6 "SSW" +1 7 "SSE" +1 8 "ESE" +2 9 "Base" +2 10 "N Even" +2 11 "N Odd" +2 12 "W Even" +2 13 "W Odd" +2 14 "S Even" +2 15 "S Odd" +2 16 "E Even" +2 17 "E Odd" +$EndPhysicalNames +$Nodes +13 +1 0 0 0 +2 0.1414213562373095 0.1414213562373095 0 +3 -0.1414213562373095 0.1414213562373095 0 +4 -0.1414213562373095 -0.1414213562373095 0 +5 0.1414213562373095 -0.1414213562373095 0 +6 1 0 0 +7 0.7071067811865475 0.7071067811865475 0 +8 0 1 0 +9 -0.7071067811865475 0.7071067811865475 0 +10 -1 0 0 +11 -0.7071067811865475 -0.7071067811865475 0 +12 0 -1 0 +13 0.7071067811865475 -0.7071067811865475 0 +$EndNodes +$Elements +20 +1 1 2 1 17 6 7 +2 1 2 2 18 7 8 +3 1 2 3 19 8 9 +4 1 2 4 20 9 10 +5 1 2 5 21 10 11 +6 1 2 6 22 11 12 +7 1 2 7 23 12 13 +8 1 2 8 24 13 6 +9 2 2 10 5 1 2 8 +10 2 2 11 6 1 8 3 +11 2 2 12 7 1 3 10 +12 2 2 13 8 1 10 4 +13 2 2 14 9 1 4 12 +14 2 2 15 10 1 12 5 +15 2 2 16 11 1 5 6 +16 2 2 17 12 1 6 2 +17 3 2 9 1 2 6 7 8 +18 3 2 9 2 3 8 9 10 +19 3 2 9 3 4 10 11 12 +20 3 2 9 4 5 12 13 6 +$EndElements diff --git a/examples/ex39.py b/examples/ex39.py new file mode 100644 index 00000000..028ff47c --- /dev/null +++ b/examples/ex39.py @@ -0,0 +1,286 @@ +''' + MFEM example 39 (converted from ex39.cpp) + + See c++ version in the MFEM library for more detail + + Sample runs: python ex39.py + python ex39.py -ess "Southern Boundary" + python ex39.py -src Base + + Description: This example code demonstrates the use of named attribute + sets in MFEM to specify material regions, boundary regions, + or source regions by name rather than attribute numbers. It + also demonstrates how new named attribute sets may be created + from arbitrary groupings of attribute numbers and used as a + convenient shorthand to refer to those groupings in other + portions of the application or through the command line. + + The particular problem being solved here is nearly the same + as that in example 1 i.e. a simple finite element + discretization of the Laplace problem -Delta u = 1 with + homogeneous Dirichlet boundary conditions and, in this case, + an inhomogeneous diffusion coefficient. The diffusion + coefficient is given a small default value throughout the + domain which is increased by two separate amounts in two named + regions. + + This example makes use of a specific input mesh, "compass.msh", + containing named domain and boundary regions generated by Gmsh + and stored in their "msh" format (version 2.2). This file + defines eight boundary regions corresponding to eight compass + headings; "ENE", "NNE", "NNW", "WSW", "SSW", "SSE", and "ESE". + It also defines nine domain regions; "Base", "N Even", "N Odd", + "W Even", "W Odd", "S Even", "S Odd", "E Even", and "E Odd". + These regions split the four compass pointers into two halves + each and also label the remaining elements as "Base". Starting + with these named regions we test the construction of named + sets as well as reading and writing these named groupings from + and to mesh files. + + The example highlights the use of named attribute sets for + both subdomains and boundaries in different contexts as well + as basic methods to create named sets from existing attributes. + +''' +import os +from os.path import expanduser, join +import numpy as np + +import mfem.ser as mfem + + +def run(order=1, + meshfile='', + source_name='Rose Even', + ess_name='Boundary', + visualization=True): + + # 2. Read the mesh from the given mesh file. We can handle triangular, + # quadrilateral, tetrahedral, hexahedral, surface and volume meshes with + # the same code. + + mesh = mfem.Mesh(meshfile, 1, 1) + dim = mesh.Dimension() + + # 3. Refine the mesh to increase the resolution. In this example we do + # 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the + # largest number that gives a final mesh with no more than 50,000 + # elements. + ref_levels = int(np.log(50000./mesh.GetNE())/np.log(2.)/dim) + for i in range(ref_levels): + mesh.UniformRefinement() + + # 4a. Display attribute set names contained in the initial mesh + # GetAttributeSetNames returns Python set object. + attr_sets = mesh.attribute_sets + bdr_attr_sets = mesh.bdr_attribute_sets + names = attr_sets.GetAttributeSetNames() + print("Element Attribute Set Names: " + + ' '.join(['"' + x + '"' for x in names])) + names = bdr_attr_sets.GetAttributeSetNames() + print("Boundary Attribute Set Names: " + + ' '.join(['"' + x + '"' for x in names])) + + # 4b. Define new regions based on existing attribute sets + Na = attr_sets.GetAttributeSet("N Even") + Nb = attr_sets.GetAttributeSet("N Odd") + Sa = attr_sets.GetAttributeSet("S Even") + Sb = attr_sets.GetAttributeSet("S Odd") + Ea = attr_sets.GetAttributeSet("E Even") + Eb = attr_sets.GetAttributeSet("E Odd") + Wa = attr_sets.GetAttributeSet("W Even") + Wb = attr_sets.GetAttributeSet("W Odd") + + # Create a new set spanning the North point + attr_sets.SetAttributeSet("North", Na) + attr_sets.AddToAttributeSet("North", Nb) + + # Create a new set spanning the South point + attr_sets.SetAttributeSet("South", Sa) + attr_sets.AddToAttributeSet("South", Sb) + + # Create a new set spanning the East point + attr_sets.SetAttributeSet("East", Ea) + attr_sets.AddToAttributeSet("East", Eb) + + # Create a new set spanning the West point + attr_sets.SetAttributeSet("West", Wa) + attr_sets.AddToAttributeSet("West", Wb) + + # Create a new set consisting of the "a" sides of the compass rose + attr_sets.SetAttributeSet("Rose Even", Na) + attr_sets.AddToAttributeSet("Rose Even", Sa) + attr_sets.AddToAttributeSet("Rose Even", Ea) + attr_sets.AddToAttributeSet("Rose Even", Wa) + + # Create a new set consisting of the "b" sides of the compass rose + attr_sets.SetAttributeSet("Rose Odd", Nb) + attr_sets.AddToAttributeSet("Rose Odd", Sb) + attr_sets.AddToAttributeSet("Rose Odd", Eb) + attr_sets.AddToAttributeSet("Rose Odd", Wb) + + # Create a new set consisting of the full compass rose + Ra = attr_sets.GetAttributeSet("Rose Even") + Rb = attr_sets.GetAttributeSet("Rose Odd") + attr_sets.SetAttributeSet("Rose", Ra) + attr_sets.AddToAttributeSet("Rose", Rb) + + # 4c. Define new boundary regions based on existing boundary attribute sets + NNE = bdr_attr_sets.GetAttributeSet("NNE") + NNW = bdr_attr_sets.GetAttributeSet("NNW") + ENE = bdr_attr_sets.GetAttributeSet("ENE") + ESE = bdr_attr_sets.GetAttributeSet("ESE") + SSE = bdr_attr_sets.GetAttributeSet("SSE") + SSW = bdr_attr_sets.GetAttributeSet("SSW") + WNW = bdr_attr_sets.GetAttributeSet("WNW") + WSW = bdr_attr_sets.GetAttributeSet("WSW") + + bdr_attr_sets.SetAttributeSet("Northern Boundary", NNE) + bdr_attr_sets.AddToAttributeSet("Northern Boundary", NNW) + + bdr_attr_sets.SetAttributeSet("Southern Boundary", SSE) + bdr_attr_sets.AddToAttributeSet("Southern Boundary", SSW) + + bdr_attr_sets.SetAttributeSet("Eastern Boundary", ENE) + bdr_attr_sets.AddToAttributeSet("Eastern Boundary", ESE) + + bdr_attr_sets.SetAttributeSet("Western Boundary", WNW) + bdr_attr_sets.AddToAttributeSet("Western Boundary", WSW) + + bdr_attr_sets.SetAttributeSet("Boundary", + bdr_attr_sets.GetAttributeSet + ("Northern Boundary")) + bdr_attr_sets.AddToAttributeSet("Boundary", + bdr_attr_sets.GetAttributeSet + ("Southern Boundary")) + bdr_attr_sets.AddToAttributeSet("Boundary", + bdr_attr_sets.GetAttributeSet + ("Eastern Boundary")) + bdr_attr_sets.AddToAttributeSet("Boundary", + bdr_attr_sets.GetAttributeSet + ("Western Boundary")) + + # 5. Define a finite element space on the mesh. Here we use continuous + # Lagrange finite elements of the specified order. + + fec = mfem.H1_FECollection(order, dim) + fespace = mfem.FiniteElementSpace(mesh, fec) + print('Number of finite element unknowns: ' + + str(fespace.GetTrueVSize())) + + # 6. Determine the list of true (i.e. conforming) essential boundary dofs. + # In this example, the boundary conditions are defined by marking all + # the boundary regions corresponding to the boundary attributes + # contained in the set named "ess_name" as essential (Dirichlet) and + # converting them to a list of true dofs. + ess_tdof_list = mfem.intArray() + if bdr_attr_sets.AttributeSetExists(ess_name): + ess_bdr_marker = bdr_attr_sets.GetAttributeSetMarker(ess_name) + fespace.GetEssentialTrueDofs(ess_bdr_marker, ess_tdof_list) + + # 7. Set up the linear form b(.) which corresponds to the right-hand side of + # the FEM linear system, which in this case is (1_s,phi_i) where phi_i + # are the basis functions in fespace and 1_s is an indicator function + # equal to 1 on the region defined by the named set "source_name" and + # zero elsewhere. + + source_marker = attr_sets.GetAttributeSetMarker(source_name) + b = mfem.LinearForm(fespace) + one = mfem.ConstantCoefficient(1.0) + b.AddDomainIntegrator(mfem.DomainLFIntegrator(one), source_marker) + b.Assemble() + + # 8. Define the solution vector x as a finite element grid function + # corresponding to fespace. Initialize x with initial guess of zero, + # which satisfies the boundary conditions. + x = mfem.GridFunction(fespace) + x.Assign(0.0) + + # 9. Set up the bilinear form a(.,.) on the finite element space + # corresponding to the Laplacian operator -Delta, by adding the + # Diffusion domain integrator. + a = mfem.BilinearForm(fespace) + + defaultCoef = mfem.ConstantCoefficient(1.0e-6) + baseCoef = mfem.ConstantCoefficient(1.0) + roseCoef = mfem.ConstantCoefficient(2.0) + + base_marker = attr_sets.GetAttributeSetMarker("Base") + rose_marker = attr_sets.GetAttributeSetMarker("Rose Even") + + # Impose a very small diffusion coefficient across the entire mesh + a.AddDomainIntegrator(mfem.DiffusionIntegrator(defaultCoef)) + + # Impose an additional, stronger diffusion coefficient in select regions + a.AddDomainIntegrator(mfem.DiffusionIntegrator(baseCoef), base_marker) + a.AddDomainIntegrator(mfem.DiffusionIntegrator(roseCoef), rose_marker) + + # 10. Assemble the bilinear form and the corresponding linear system, + # applying any necessary transformations. + a.Assemble() + + A = mfem.SparseMatrix() + B = mfem.Vector() + X = mfem.Vector() + a.FormLinearSystem(ess_tdof_list, x, b, A, X, B) + print("Size of linear system: " + str(A.Height())) + + # 11. Solve the system using PCG with symmetric Gauss-Seidel preconditioner. + M = mfem.GSSmoother(A) + mfem.PCG(A, M, B, X, 1, 800, 1e-12, 0.0) + + # 12. Recover the solution as a finite element grid function. + a.RecoverFEMSolution(X, b, x) + + # 13. Save the refined mesh and the solution. This output can be viewed + # later using GLVis: "glvis -m refined.mesh -g sol.gf". + mesh.Save('refined.mesh') + x.Save('sol.gf') + + + # 14. Send the solution by socket to a GLVis server. + if visualization: + sol_sock = mfem.socketstream("localhost", 19916) + sol_sock.precision(8) + sol_sock << "solution\n" << mesh << x << "keys Rjmm" + sol_sock.flush() + + +if __name__ == "__main__": + from mfem.common.arg_parser import ArgParser + + parser = ArgParser(description='Ex1 (Laplace Problem)') + parser.add_argument('-m', '--mesh', + default='compass.msh', + action='store', type=str, + help='Mesh file to use.') + parser.add_argument('-o', '--order', + action='store', default=1, type=int, + help='Finite element order (polynomial degree) or -1 for isoparametric space.') + parser.add_argument('-src', '--source-attr-name', + action='store', default='Rose Even', type=str, + help='Name of attribute set containing source.') + parser.add_argument('-ess', '--ess-attr-name', + action='store', default='Boundary', type=str, + help='Name of attribute set containing essential BC.') + parser.add_argument('-no-vis', '--no-visualization', + action='store_true', + default=False, + help='Disable or disable GLVis visualization') + + args = parser.parse_args() + parser.print_options(args) + + order = args.order + meshfile = expanduser( + join(os.path.dirname(__file__), '..', 'data', args.mesh)) + + visualization = not args.no_visualization + source_name = args.source_attr_name + ess_name = args.ess_attr_name + + run(order=order, + meshfile=meshfile, + source_name=source_name, + ess_name=ess_name, + visualization=visualization) diff --git a/examples/ex39p.py b/examples/ex39p.py new file mode 100644 index 00000000..3756203e --- /dev/null +++ b/examples/ex39p.py @@ -0,0 +1,309 @@ +''' + MFEM example 39 (converted from ex39.cpp) + + See c++ version in the MFEM library for more detail + + Sample runs: mpirun -np 4 python ex39p.py + mpirun -np 4 python ex39p.py -ess "Southern Boundary" + mpirun -np 4 python ex39p.py -src Base + + Description: This example code demonstrates the use of named attribute + sets in MFEM to specify material regions, boundary regions, + or source regions by name rather than attribute numbers. It + also demonstrates how new named attribute sets may be created + from arbitrary groupings of attribute numbers and used as a + convenient shorthand to refer to those groupings in other + portions of the application or through the command line. + + The particular problem being solved here is nearly the same + as that in example 1 i.e. a simple finite element + discretization of the Laplace problem -Delta u = 1 with + homogeneous Dirichlet boundary conditions and, in this case, + an inhomogeneous diffusion coefficient. The diffusion + coefficient is given a small default value throughout the + domain which is increased by two separate amounts in two named + regions. + + This example makes use of a specific input mesh, "compass.msh", + containing named domain and boundary regions generated by Gmsh + and stored in their "msh" format (version 2.2). This file + defines eight boundary regions corresponding to eight compass + headings; "ENE", "NNE", "NNW", "WSW", "SSW", "SSE", and "ESE". + It also defines nine domain regions; "Base", "N Even", "N Odd", + "W Even", "W Odd", "S Even", "S Odd", "E Even", and "E Odd". + These regions split the four compass pointers into two halves + each and also label the remaining elements as "Base". Starting + with these named regions we test the construction of named + sets as well as reading and writing these named groupings from + and to mesh files. + + The example highlights the use of named attribute sets for + both subdomains and boundaries in different contexts as well + as basic methods to create named sets from existing attributes. + +''' +import os +from os.path import expanduser, join +import numpy as np + +import mfem.par as mfem + +from mpi4py import MPI +num_procs = MPI.COMM_WORLD.size +myid = MPI.COMM_WORLD.rank + + +def run(order=1, + meshfile='', + source_name='Rose Even', + ess_name='Boundary', + visualization=True): + + # 3. Read the serial mesh from the given mesh file. + mesh = mfem.Mesh(meshfile, 1, 1) + dim = mesh.Dimension() + + # 4. Refine the serial mesh on all processors to increase the resolution. In + # this example we do 'ref_levels' of uniform refinement. We choose + # 'ref_levels' to be the largest number that gives a final mesh with no + # more than 10,000 elements. + + ref_levels = int(np.log(10000./mesh.GetNE())/np.log(2.)/dim) + for i in range(ref_levels): + mesh.UniformRefinement() + + # 5. Define a parallel mesh by a partitioning of the serial mesh. Refine + # this mesh further in parallel to increase the resolution. Once the + # parallel mesh is defined, the serial mesh can be deleted. + pmesh = mfem.ParMesh(MPI.COMM_WORLD, mesh) + del mesh + + par_ref_levels = 2 + for l in range(par_ref_levels): + pmesh.UniformRefinement() + + # 6a. Display attribute set names contained in the initial mesh + # GetAttributeSetNames returns Python set object. + + attr_sets = pmesh.attribute_sets + bdr_attr_sets = pmesh.bdr_attribute_sets + + if myid == 0: + names = attr_sets.GetAttributeSetNames() + print("Element Attribute Set Names: " + + ' '.join(['"' + x + '"' for x in names])) + names = bdr_attr_sets.GetAttributeSetNames() + print("Boundary Attribute Set Names: " + + ' '.join(['"' + x + '"' for x in names])) + + # 6b. Define new regions based on existing attribute sets + Na = attr_sets.GetAttributeSet("N Even") + Nb = attr_sets.GetAttributeSet("N Odd") + Sa = attr_sets.GetAttributeSet("S Even") + Sb = attr_sets.GetAttributeSet("S Odd") + Ea = attr_sets.GetAttributeSet("E Even") + Eb = attr_sets.GetAttributeSet("E Odd") + Wa = attr_sets.GetAttributeSet("W Even") + Wb = attr_sets.GetAttributeSet("W Odd") + + # Create a new set spanning the North point + attr_sets.SetAttributeSet("North", Na) + attr_sets.AddToAttributeSet("North", Nb) + + # Create a new set spanning the South point + attr_sets.SetAttributeSet("South", Sa) + attr_sets.AddToAttributeSet("South", Sb) + + # Create a new set spanning the East point + attr_sets.SetAttributeSet("East", Ea) + attr_sets.AddToAttributeSet("East", Eb) + + # Create a new set spanning the West point + attr_sets.SetAttributeSet("West", Wa) + attr_sets.AddToAttributeSet("West", Wb) + + # Create a new set consisting of the "a" sides of the compass rose + attr_sets.SetAttributeSet("Rose Even", Na) + attr_sets.AddToAttributeSet("Rose Even", Sa) + attr_sets.AddToAttributeSet("Rose Even", Ea) + attr_sets.AddToAttributeSet("Rose Even", Wa) + + # Create a new set consisting of the "b" sides of the compass rose + attr_sets.SetAttributeSet("Rose Odd", Nb) + attr_sets.AddToAttributeSet("Rose Odd", Sb) + attr_sets.AddToAttributeSet("Rose Odd", Eb) + attr_sets.AddToAttributeSet("Rose Odd", Wb) + + # Create a new set consisting of the full compass rose + Ra = attr_sets.GetAttributeSet("Rose Even") + Rb = attr_sets.GetAttributeSet("Rose Odd") + attr_sets.SetAttributeSet("Rose", Ra) + attr_sets.AddToAttributeSet("Rose", Rb) + + # 6c. Define new boundary regions based on existing boundary attribute sets + NNE = bdr_attr_sets.GetAttributeSet("NNE") + NNW = bdr_attr_sets.GetAttributeSet("NNW") + ENE = bdr_attr_sets.GetAttributeSet("ENE") + ESE = bdr_attr_sets.GetAttributeSet("ESE") + SSE = bdr_attr_sets.GetAttributeSet("SSE") + SSW = bdr_attr_sets.GetAttributeSet("SSW") + WNW = bdr_attr_sets.GetAttributeSet("WNW") + WSW = bdr_attr_sets.GetAttributeSet("WSW") + + bdr_attr_sets.SetAttributeSet("Northern Boundary", NNE) + bdr_attr_sets.AddToAttributeSet("Northern Boundary", NNW) + + bdr_attr_sets.SetAttributeSet("Southern Boundary", SSE) + bdr_attr_sets.AddToAttributeSet("Southern Boundary", SSW) + + bdr_attr_sets.SetAttributeSet("Eastern Boundary", ENE) + bdr_attr_sets.AddToAttributeSet("Eastern Boundary", ESE) + + bdr_attr_sets.SetAttributeSet("Western Boundary", WNW) + bdr_attr_sets.AddToAttributeSet("Western Boundary", WSW) + + bdr_attr_sets.SetAttributeSet("Boundary", + bdr_attr_sets.GetAttributeSet + ("Northern Boundary")) + bdr_attr_sets.AddToAttributeSet("Boundary", + bdr_attr_sets.GetAttributeSet + ("Southern Boundary")) + bdr_attr_sets.AddToAttributeSet("Boundary", + bdr_attr_sets.GetAttributeSet + ("Eastern Boundary")) + bdr_attr_sets.AddToAttributeSet("Boundary", + bdr_attr_sets.GetAttributeSet + ("Western Boundary")) + + # 7. Define a parallel finite element space on the parallel mesh. Here we + # use continuous Lagrange finite elements of the specified order. If + # order < 1, we instead use an isoparametric/isogeometric space. + fec = mfem.H1_FECollection(order, dim) + fespace = mfem.ParFiniteElementSpace(pmesh, fec) + size = fespace.GlobalTrueVSize() + if myid == 0: + print('Number of finite element unknowns: ' + str(size)) + + # 8. Determine the list of true (i.e. parallel conforming) essential + # boundary dofs. In this example, the boundary conditions are defined + # by marking all the boundary regions corresponding to the boundary + # attributes contained in the set named "ess_name" as essential + # (Dirichlet) and converting them to a list of true dofs. + ess_tdof_list = mfem.intArray() + if bdr_attr_sets.AttributeSetExists(ess_name): + ess_bdr_marker = bdr_attr_sets.GetAttributeSetMarker(ess_name) + fespace.GetEssentialTrueDofs(ess_bdr_marker, ess_tdof_list) + + # 9. Set up the parallel linear form b(.) which corresponds to the + # right-hand side of the FEM linear system, which in this case is + # (1_s,phi_i) where phi_i are the basis functions in fespace and 1_s + # is an indicator function equal to 1 on the region defined by the + # named set "source_name" and zero elsewhere. + + source_marker = attr_sets.GetAttributeSetMarker(source_name) + b = mfem.ParLinearForm(fespace) + one = mfem.ConstantCoefficient(1.0) + b.AddDomainIntegrator(mfem.DomainLFIntegrator(one), source_marker) + b.Assemble() + + # 10. Define the solution vector x as a parallel finite element grid + # function corresponding to fespace. Initialize x with initial guess of + # zero, which satisfies the boundary conditions. + x = mfem.ParGridFunction(fespace) + x.Assign(0.0) + + # 11. Set up the parallel bilinear form a(.,.) on the finite element space + # corresponding to the Laplacian operator -Delta, by adding the + # Diffusion domain integrator. + a = mfem.ParBilinearForm(fespace) + + defaultCoef = mfem.ConstantCoefficient(1.0e-6) + baseCoef = mfem.ConstantCoefficient(1.0) + roseCoef = mfem.ConstantCoefficient(2.0) + + base_marker = attr_sets.GetAttributeSetMarker("Base") + rose_marker = attr_sets.GetAttributeSetMarker("Rose Even") + + # Impose a very small diffusion coefficient across the entire mesh + a.AddDomainIntegrator(mfem.DiffusionIntegrator(defaultCoef)) + + # Impose an additional, stronger diffusion coefficient in select regions + a.AddDomainIntegrator(mfem.DiffusionIntegrator(baseCoef), base_marker) + a.AddDomainIntegrator(mfem.DiffusionIntegrator(roseCoef), rose_marker) + + # 12. Assemble the bilinear form and the corresponding linear system, + # applying any necessary transformations. + a.Assemble() + + A = mfem.HypreParMatrix() + B = mfem.Vector() + X = mfem.Vector() + a.FormLinearSystem(ess_tdof_list, x, b, A, X, B) + + # 13. Solve the system using PCG with hypre's BoomerAMG preconditioner. + M = mfem.HypreBoomerAMG(A) + cg = mfem.CGSolver(MPI.COMM_WORLD) + cg.SetRelTol(1e-12) + cg.SetMaxIter(2000) + cg.SetPrintLevel(1) + cg.SetPreconditioner(M) + cg.SetOperator(A) + cg.Mult(B, X) + + # 14. Recover the parallel grid function corresponding to X. This is the + # local finite element solution on each processor. + a.RecoverFEMSolution(X, b, x) + + # 15. Save the refined mesh and the solution in parallel. This output can + # be viewed later using GLVis: "glvis -np -m mesh -g sol". + pmesh.Save('mesh') + x.Save('sol') + + # 16. Send the solution by socket to a GLVis server. + if visualization: + sol_sock = mfem.socketstream("localhost", 19916) + sol_sock << "parallel " << num_procs << " " << myid << "\n" + sol_sock.precision(8) + sol_sock << "solution\n" << pmesh << x << "keys Rjmm" + sol_sock.flush() + + +if __name__ == "__main__": + from mfem.common.arg_parser import ArgParser + + parser = ArgParser(description='Ex1 (Laplace Problem)') + parser.add_argument('-m', '--mesh', + default='compass.msh', + action='store', type=str, + help='Mesh file to use.') + parser.add_argument('-o', '--order', + action='store', default=1, type=int, + help='Finite element order (polynomial degree) or -1 for isoparametric space.') + parser.add_argument('-src', '--source-attr-name', + action='store', default='Rose Even', type=str, + help='Name of attribute set containing source.') + parser.add_argument('-ess', '--ess-attr-name', + action='store', default='Boundary', type=str, + help='Name of attribute set containing essential BC.') + parser.add_argument('-no-vis', '--no-visualization', + action='store_true', + default=False, + help='Disable or disable GLVis visualization') + + args = parser.parse_args() + if myid == 0: + parser.print_options(args) + + order = args.order + meshfile = expanduser( + join(os.path.dirname(__file__), '..', 'data', args.mesh)) + + visualization = not args.no_visualization + source_name = args.source_attr_name + ess_name = args.ess_attr_name + + run(order=order, + meshfile=meshfile, + source_name=source_name, + ess_name=ess_name, + visualization=visualization) diff --git a/mfem/_par/attribute_sets.i b/mfem/_par/attribute_sets.i index c151dc3a..1f130ccf 100644 --- a/mfem/_par/attribute_sets.i +++ b/mfem/_par/attribute_sets.i @@ -24,6 +24,18 @@ import_array(); OSTREAM_TYPEMAP(std::ostream&) ISTREAM_TYPEMAP(std::istream&) +// +// AttributeSets::GetAttributeSetNames returns Python set +// +%typemap(out) std::set{ + resultobj = PySet_New(NULL); + for (auto const &set_name : *(& $1)){ + std::ostringstream name2; + name2 << set_name; + PySet_Add(resultobj, PyUnicode_FromString(name2.str().c_str())); + } +} + %include "mesh/attribute_sets.hpp" #endif diff --git a/mfem/_ser/attribute_sets.i b/mfem/_ser/attribute_sets.i index 6729597b..b9fd738f 100644 --- a/mfem/_ser/attribute_sets.i +++ b/mfem/_ser/attribute_sets.i @@ -1,4 +1,6 @@ %module(package="mfem._ser") attribute_sets +%feature("autodoc", "1"); + %{ #include "mfem.hpp" #include "numpy/arrayobject.h" @@ -16,6 +18,7 @@ import_array(); %} %include "exception.i" +%include "std_string.i" %include "../common/exception.i" %import "array.i" @@ -24,6 +27,21 @@ import_array(); OSTREAM_TYPEMAP(std::ostream&) ISTREAM_TYPEMAP(std::istream&) +// +// AttributeSets::GetAttributeSetNames returns Python set +// +%typemap(out) std::set{ + resultobj = PySet_New(NULL); + for (auto const &set_name : *(& $1)){ + std::ostringstream name2; + name2 << set_name; + PySet_Add(resultobj, PyUnicode_FromString(name2.str().c_str())); + } +} + %include "mesh/attribute_sets.hpp" #endif + + + diff --git a/mfem/par.py b/mfem/par.py index 47b667c7..bcf83e1d 100644 --- a/mfem/par.py +++ b/mfem/par.py @@ -77,6 +77,7 @@ from mfem._par.qfunction import * from mfem._par.quadinterpolator import * from mfem._par.quadinterpolator_face import * +from mfem._par.attribute_sets import * from mfem._par.fe_base import * from mfem._par.fe_h1 import * diff --git a/mfem/ser.py b/mfem/ser.py index 3ff0d3ca..ba28e29e 100644 --- a/mfem/ser.py +++ b/mfem/ser.py @@ -60,6 +60,7 @@ from mfem._ser.qfunction import * from mfem._ser.quadinterpolator import * from mfem._ser.quadinterpolator_face import * +from mfem._ser.attribute_sets import * from mfem._ser.fe_base import * from mfem._ser.fe_h1 import *