@@ -49,9 +49,8 @@ fn main() -> anyhow::Result<()> {
4949 // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
5050 // largest number that gives a final mesh with no more than 50,000
5151 // elements.
52- let ref_levels = ( ( 50000. / mesh. get_num_elems ( ) as f64 ) . log2 ( )
53- / dim as f64 )
54- . floor ( ) as u32 ;
52+ let ne = mesh. get_num_elems ( ) as f64 ;
53+ let ref_levels = ( ( 50000. / ne) . log2 ( ) / dim as f64 ) . floor ( ) as u32 ;
5554 for _ in 0 ..ref_levels {
5655 mesh. uniform_refinement ( RefAlgo :: A ) ;
5756 }
@@ -60,33 +59,17 @@ fn main() -> anyhow::Result<()> {
6059 // 5. Define a finite element space on the mesh. Here we use continuous
6160 // Lagrange finite elements of the specified order. If order < 1, we
6261 // instead use an isoparametric/isogeometric space.
63- let owned_fec: Option < H1FeCollection > = if args. order > 0 {
64- Some ( H1FeCollection :: new (
65- args. order ,
66- dim,
67- BasisType :: GaussLobatto ,
68- ) )
69- } else if mesh. get_nodes ( ) . is_none ( ) {
70- Some ( H1FeCollection :: new ( 1 , dim, BasisType :: GaussLobatto ) )
71- } else {
72- None
73- } ;
74-
75- let owned_nodes = mesh. get_nodes ( ) ;
76-
77- let fec: & dyn FiniteElementCollection = match & owned_fec {
78- Some ( h1_fec) => h1_fec,
79- None => {
80- println ! ( "Using isoparametric FEs" ) ;
81- let nodes = owned_nodes. as_ref ( ) . expect ( "Mesh has its own nodes" ) ;
82- let iso_fec = nodes. get_own_fec ( ) . expect ( "OwnFEC exists" ) ;
83- iso_fec
62+ let mut mesh_fec = mesh. with_fec ( |fec| {
63+ if args. order > 0 {
64+ H1_FECollection :: new ( args. order , dim) . into ( )
65+ } else if let Some ( fec) = fec {
66+ fec. into ( )
67+ } else {
68+ H1_FECollection :: new ( 1 , dim) . into ( )
8469 }
85- } ;
86-
87- dbg ! ( fec. get_name( ) ) ;
88-
89- let fespace = FiniteElementSpace :: new ( & mesh, fec, 1 , OrderingType :: byNODES) ;
70+ } ) ;
71+ dbg ! ( mesh_fec. fec( ) . get_name( ) ) ;
72+ let fespace = FiniteElementSpace :: new ( & mut mesh_fec) . build ( ) ;
9073 println ! (
9174 "Number of finite element unknowns: {}" ,
9275 fespace. get_true_vsize( ) ,
@@ -97,33 +80,32 @@ fn main() -> anyhow::Result<()> {
9780 // the boundary attributes from the mesh as essential (Dirichlet) and
9881 // converting them to a list of true dofs.
9982 let mut ess_tdof_list = ArrayInt :: new ( ) ;
100- if let Some ( max_bdr_attr) = mesh. get_bdr_attributes ( ) . iter ( ) . max ( ) {
83+ if let Some ( max_bdr_attr) = fespace . mesh ( ) . bdr_attributes ( ) . iter ( ) . max ( ) {
10184 let mut ess_bdr = ArrayInt :: with_len ( * max_bdr_attr as usize ) ;
102- ess_bdr. set_all ( 1 ) ;
85+ ess_bdr. fill ( 1 ) ;
10386 fespace. get_essential_true_dofs ( & ess_bdr, & mut ess_tdof_list, None ) ;
10487 }
10588
106- // 7. Set up the linear form b(.) which corresponds to the right-hand side of
107- // the FEM linear system, which in this case is (1,phi_i) where phi_i are
108- // the basis functions in the finite element fespace.
89+ // 7. Set up the linear form b(.) which corresponds to the
90+ // right-hand side of the FEM linear system, which in this case
91+ // is (1,phi_i) where phi_i are the basis functions in the
92+ // finite element `fespace`.
93+ let mut one = ConstantCoefficient :: new ( 1.0 ) ;
10994 let mut b = LinearForm :: new ( & fespace) ;
110- let one = ConstantCoefficient :: new ( 1.0 ) ;
111- let integrator = DomainLFIntegrator :: new ( & one, 2 , 0 ) ;
112- b. add_domain_integrator ( integrator) ;
95+ b. add_domain_integrator ( DomainLFIntegrator :: new ( & mut one) ) ;
11396 b. assemble ( ) ;
11497
11598 // 8. Define the solution vector x as a finite element grid function
11699 // corresponding to fespace. Initialize x with initial guess of zero,
117100 // which satisfies the boundary conditions.
118101 let mut x = GridFunction :: new ( & fespace) ;
119- x. set_all ( 0.0 ) ;
102+ x. fill ( 0.0 ) ;
120103
121104 // 9. Set up the bilinear form a(.,.) on the finite element space
122- // corresponding to the Laplacian operator -Delta, by adding the Diffusion
123- // domain integrator.
105+ // corresponding to the Laplacian operator -Delta, by adding
106+ // the Diffusion domain integrator.
124107 let mut a = BilinearForm :: new ( & fespace) ;
125- let bf_integrator = DiffusionIntegrator :: new ( & one) ;
126- a. add_domain_integrator ( bf_integrator) ;
108+ a. add_domain_integrator ( DiffusionIntegrator :: new ( ) ) ;
127109
128110 // 10. Assemble the bilinear form and the corresponding linear system,
129111 // applying any necessary transformations such as: eliminating boundary
@@ -136,30 +118,42 @@ fn main() -> anyhow::Result<()> {
136118 let mut x_vec = Vector :: new ( ) ;
137119 a. form_linear_system (
138120 & ess_tdof_list,
139- & x,
140- & b,
121+ & mut x,
122+ & mut b,
141123 & mut a_mat,
142124 & mut x_vec,
143125 & mut b_vec,
126+ false ,
144127 ) ;
145-
146- println ! ( "Size of linear system: {}" , a_mat. height( ) ) ;
128+ println ! ( "Size of linear system: {}" , a_mat. op( ) . height( ) ) ;
147129 dbg ! ( a_mat. get_type( ) ) ;
148130
149131 // 11. Solve the linear system A X = B.
150- // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
151- let a_sparse =
152- SparseMatrixRef :: try_from ( & a_mat) . expect ( "Operator is a SparseMatrix" ) ;
153- let mut m_mat = GsSmoother :: new ( & a_sparse, 0 , 1 ) ;
154- solve_with_pcg ( & a_mat, & mut m_mat, & b_vec, & mut x_vec, 1 , 200 , 1e-12 , 0.0 ) ;
132+ // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
133+ let a_sparse: Ref < ' _ , SparseMatrix > = ( & * a_mat) . try_into ( ) ?;
134+ let mut m_mat = GSSmoother :: with_matrix ( & a_sparse, 0 , 1 ) ;
135+ mfem:: pcg (
136+ & a_mat. op ( ) ,
137+ & mut m_mat,
138+ & b_vec,
139+ & mut x_vec,
140+ 1 ,
141+ 200 ,
142+ 1e-12 ,
143+ 0.0 ,
144+ ) ;
155145
156146 // 12. Recover the solution as a finite element grid function.
157147 a. recover_fem_solution ( & x_vec, & b, & mut x) ;
158148
159- // 13. Save the refined mesh and the solution. This output can be viewed later
160- // using GLVis: "glvis -m refined.mesh -g sol.gf".
161- mesh. save_to_file ( "refined.mesh" , 8 ) ;
162- x. save_to_file ( "sol.gf" , 8 ) ;
149+ // 13. Save the refined mesh and the solution. This output can be
150+ // viewed later using GLVis: "glvis -m refined.mesh -g sol.gf".
151+ fespace . mesh ( ) . save ( ) . to_file ( "refined.mesh" ) ;
152+ x. save ( ) . to_file ( "sol.gf" ) ;
163153
164154 Ok ( ( ) )
165155}
156+
157+ // Local Variables:
158+ // compile-command: "cargo run --example ex1 -- --mesh data/square.mesh"
159+ // End:
0 commit comments