Second-order Quadrilaterals (2D) meshes in standard Abaqus .inp for…

…mat (#2217) * First steps towards quadratic/curved elements form standard Abaqus * idea * Notes * preproc abaqus * shorten * work on quad mesches * Seemingly working version 2nd order Abaqus 2D * Comments * experiment with viz * test hohqmesh * revert * comments * SD7003 * start 3d * continue * Continue * Reduce to 2D * example * todo * Mention that boundary is in general only higher-order, not curved * rename * test + example * docs * Comment * reset project toml * news * typo * Update src/meshes/p4est_mesh.jl * Update NEWS.md * Update src/meshes/p4est_mesh.jl Co-authored-by: Andrew Winters <[email protected]> * Update examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/meshes/p4est_mesh.jl * comments * shorten * try to make MPI rdy * work on mpi * datatype * test hybrid mesh --------- Co-authored-by: Andrew Winters <[email protected]> Co-authored-by: Hendrik Ranocha <[email protected]>
trixi-framework · Jan 12, 2025 · e93de2a · e93de2a
1 parent 7b7fa32
commit e93de2a
Show file tree

Hide file tree

Showing 9 changed files with 721 additions and 51 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -13,6 +13,7 @@ for human readability.
   and [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) ([#2008])
 - `LobattoLegendreBasis` and related datastructures made fully floating-type general,
   enabling calculations with higher than double (`Float64`) precision ([#2128])
+- In 2D, quadratic elements, i.e., 8-node (quadratic) quadrilaterals are now supported in standard Abaqus `inp` format ([#2217])
 
 #### Changed
 

diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl
@@ -315,7 +315,7 @@ end #hide #md
 # ```
 # which forces `gmsh` to generate quadrilateral elements instead of the default triangles.
 # This is strictly required to be able to use the mesh later with `p4est`, which supports only straight-sided quads,
-# i.e., `C2D4, CPS4, S4` in 2D and `C3D` in 3D.
+# i.e., `CPS4, S4, ...` in 2D and `C3D8` in 3D.
 # See for more details the (short) [documentation](https://p4est.github.io/p4est-howto.pdf) on the interaction of `p4est` with `.inp` files.
 # In principle, it should also be possible to use the `recombine` function of `gmsh` to convert the triangles to quads, 
 # but this is observed to be less robust than enforcing quads from the beginning.

diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md
@@ -386,6 +386,7 @@ transfinite map of the straight sided hexahedral element to find
 
 Also for a mesh in standard Abaqus format there are no qualitative changes when going from 2D to 3D.
 The most notable difference is that boundaries are formed in 3D by faces defined by four nodes while in 2D boundaries are edges consisting of two elements.
+Note that standard Abaqus also defines quadratic element types. In Trixi.jl, these higher-order elements are currently only supported in 2D, i.e., 8-node quadrilaterals.
 A simple mesh file, which is used also in `examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl`, is given below:
 ```
 *Heading

diff --git a/examples/p4est_2d_dgsem/elixir_euler_free_stream_hybrid_mesh.jl b/examples/p4est_2d_dgsem/elixir_euler_free_stream_hybrid_mesh.jl
@@ -0,0 +1,59 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+# Test free stream preservation with constant initial condition
+initial_condition = initial_condition_constant
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+###############################################################################
+
+# Hybrid mesh composed of a second-order and first-order quadrilateral element
+mesh_file = Trixi.download("https://gist.githubusercontent.com/DanielDoehring/84d876776e379322c38e9c0bfcadee8e/raw/0068805b853a8faa0fe229280d353016359d8a7d/hybrid_quadmesh.inp",
+                           joinpath(@__DIR__, "hybrid_quadmesh.inp"))
+
+# Refine the given mesh twice
+mesh = P4estMesh{2}(mesh_file, initial_refinement_level = 2)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions =
+                                    Dict(:all => BoundaryConditionDirichlet(initial_condition)))
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 2.0)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl
@@ -0,0 +1,151 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+gamma = 1.4
+
+U_inf = 0.2
+aoa = 4 * pi / 180
+rho_inf = 1.4 # with gamma = 1.4 => p_inf = 1.0
+
+Re = 10000.0
+airfoil_cord_length = 1.0
+
+t_c = airfoil_cord_length / U_inf
+
+prandtl_number() = 0.72
+mu() = rho_inf * U_inf * airfoil_cord_length / Re
+
+equations = CompressibleEulerEquations2D(gamma)
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(),
+                                                          Prandtl = prandtl_number(),
+                                                          gradient_variables = GradientVariablesPrimitive())
+
+@inline function initial_condition_mach02_flow(x, t, equations)
+    # set the freestream flow parameters such that c_inf = 1.0 => Mach 0.2
+    rho_freestream = 1.4
+
+    # Values correspond to `aoa = 4 * pi / 180`
+    v1 = 0.19951281005196486 # 0.2 * cos(aoa)
+    v2 = 0.01395129474882506 # 0.2 * sin(aoa)
+
+    p_freestream = 1.0
+
+    prim = SVector(rho_freestream, v1, v2, p_freestream)
+    return prim2cons(prim, equations)
+end
+initial_condition = initial_condition_mach02_flow
+
+surf_flux = flux_hllc
+vol_flux = flux_chandrashekar
+solver = DGSEM(polydeg = 3, surface_flux = surf_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(vol_flux))
+
+###############################################################################
+# Get the uncurved mesh from a file (downloads the file if not available locally)
+
+# Get quadratic meshfile
+mesh_file_name = "SD7003_2D_Quadratic.inp"
+
+mesh_file = Trixi.download("https://gist.githubusercontent.com/DanielDoehring/bd2aa20f7e6839848476a0e87ede9f69/raw/1bc8078b4a57634819dc27010f716e68a225c9c6/SD7003_2D_Quadratic.inp",
+                           joinpath(@__DIR__, mesh_file_name))
+
+# There is also a linear mesh file available at
+# https://gist.githubusercontent.com/DanielDoehring/375df933da8a2081f58588529bed21f0/raw/18592aa90f1c86287b4f742fd405baf55c3cf133/SD7003_2D_Linear.inp
+
+boundary_symbols = [:Airfoil, :FarField]
+mesh = P4estMesh{2}(mesh_file, boundary_symbols = boundary_symbols)
+
+boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition)
+
+velocity_bc_airfoil = NoSlip((x, t, equations) -> SVector(0.0, 0.0))
+heat_bc = Adiabatic((x, t, equations) -> 0.0)
+boundary_condition_airfoil = BoundaryConditionNavierStokesWall(velocity_bc_airfoil, heat_bc)
+
+boundary_conditions_hyp = Dict(:FarField => boundary_condition_free_stream,
+                               :Airfoil => boundary_condition_slip_wall)
+
+boundary_conditions_para = Dict(:FarField => boundary_condition_free_stream,
+                                :Airfoil => boundary_condition_airfoil)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions = (boundary_conditions_hyp,
+                                                                    boundary_conditions_para))
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Run simulation until initial pressure wave is gone.
+# Note: This is a very long simulation!
+tspan = (0.0, 30 * t_c)
+
+# Drag/Lift coefficient measurements should then be done over the 30 to 35 t_c interval
+# by restarting the simulation.
+
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+f_aoa() = aoa
+f_rho_inf() = rho_inf
+f_U_inf() = U_inf
+f_linf() = airfoil_cord_length
+
+drag_coefficient = AnalysisSurfaceIntegral((:Airfoil,),
+                                           DragCoefficientPressure(f_aoa(), f_rho_inf(),
+                                                                   f_U_inf(), f_linf()))
+
+drag_coefficient_shear_force = AnalysisSurfaceIntegral((:Airfoil,),
+                                                       DragCoefficientShearStress(f_aoa(),
+                                                                                  f_rho_inf(),
+                                                                                  f_U_inf(),
+                                                                                  f_linf()))
+
+lift_coefficient = AnalysisSurfaceIntegral((:Airfoil,),
+                                           LiftCoefficientPressure(f_aoa(), f_rho_inf(),
+                                                                   f_U_inf(), f_linf()))
+
+# For long simulation run, use a large interval.
+# For measurements once the simulation has settled in, one should use a 
+# significantly smaller interval, e.g. 500 to record the drag/lift coefficients.                                                                   
+analysis_interval = 10_000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     output_directory = "out",
+                                     save_analysis = true,
+                                     analysis_errors = Symbol[], # Turn off standard errors
+                                     analysis_integrals = (drag_coefficient,
+                                                           drag_coefficient_shear_force,
+                                                           lift_coefficient))
+
+stepsize_callback = StepsizeCallback(cfl = 2.2)
+
+alive_callback = AliveCallback(alive_interval = 50)
+
+save_solution = SaveSolutionCallback(interval = analysis_interval,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim,
+                                     output_directory = "out")
+
+save_restart = SaveRestartCallback(interval = analysis_interval,
+                                   save_final_restart = true)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        stepsize_callback,
+                        save_solution,
+                        save_restart);
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode,
+            CarpenterKennedy2N54(williamson_condition = false,
+                                 thread = OrdinaryDiffEq.True());
+            dt = 1.0, save_everystep = false, callback = callbacks)
+
+summary_callback() # print the timer summary