diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index abcfcf8f982..1e5d782ebb6 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -3,7 +3,7 @@ The [`P4estMesh`](@ref) is an unstructured, curvilinear, nonconforming mesh type for quadrilateral (2D) and hexahedral (3D) cells. It supports quadtree/octree-based adaptive mesh refinement (AMR) via -the C library [p4est](https://github.com/cburstedde/p4est). See +the C library [`p4est`](https://github.com/cburstedde/p4est). See [`AMRCallback`](@ref) for further information. Due to its curvilinear nature, (numerical) fluxes need to implement methods @@ -12,3 +12,338 @@ such as the compressible Euler equations can use [`FluxRotated`](@ref) to wrap numerical fluxes implemented only for Cartesian meshes. This simplifies the re-use of existing functionality for the [`TreeMesh`](@ref) but is usually less efficient, cf. [PR #550](https://github.com/trixi-framework/Trixi.jl/pull/550). + +## Construction of a P4estMesh from an Abaqus file + +One available option to construct a [`P4estMesh`](@ref) is to read in an Abaqus (`.inp`) mesh file. +We briefly describe the structure of this file, the conventions it uses, and how the mesh file +is parsed to create an initial unstructured, curvilinear, and conforming mesh. + +### Mesh in two spatial dimensions + +For this discussion we use the following two-dimensional unstructured curved mesh with three elements: + +![abaqus-2dmesh-docs](https://user-images.githubusercontent.com/25242486/139241997-88e70a01-286f-4cee-80b1-2fd83c60bcca.png) + +We note that the corner and element connectivity information parsed from the Abaqus file creates +a straight sided (linear) mesh. +From this linear mesh there are two strategies available to make the mesh curvilinear: + +1. Apply a `mapping` function to describe a transformation of the linear mesh to another + physical domain. The mapping is approximated using interpolation polynomial of a user + specified polynomial degree. The default value of this polynomial degree is `1` that + corresponds to an uncurved geometry. +2. High-order boundary information is available in the `.inp` mesh file because it was created + with the [HOHQMesh](https://github.com/trixi-framework/HOHQMesh) mesh generator, which + is available via the Julia package [HOHQMesh.jl](https://github.com/trixi-framework/HOHQMesh.jl). + This information is used to create appropriate transfinite mappings during the mesh construction. + +We divide our discussion into two parts. The first part discusses the standard corner and element information +contained in the `.inp` mesh file. The second part specifically deals with the mesh file parsing of an Abaqus +file created by HOHQMesh.jl. + +#### Mesh file header + +An Abaqus `.inp` mesh file typically begins with a `*Heading`. +Though *optional*, the `*Heading` is helpful to give users some information about the mesh described by the mesh file. +In particular, a `.inp` mesh file created with `HOHQMesh` will contain the header +``` +*Heading + File created by HOHQMesh +``` +This heading is used to indicate to the mesh constructor which of the above mapping strategies to apply in order to +create a curvilinear mesh. +If the Abaqus file header is **not** present then the `P4estMesh` is created with the first strategy above. + +#### List of corner nodes + +Next, prefaced with `*NODE`, comes a list of the physical `(x,y,z)` coordinates of all the corners. +The first integer in the list of the corners provides its id number. +Thus, for the two-dimensional example mesh this block of corner information is +``` +*NODE +1, 1.0, -1.0, 0.0 +2, 3.0, 0.0, 0.0 +3, 1.0, 1.0, 0.0 +4, 2.0, 0.0, 0.0 +5, 0.0, 0.0, 0.0 +6, 3.0, 1.0, 0.0 +7, 3.0, -1.0, 0.0 +``` + +#### List of elements + +The element connectivity is given after the list of corners. The header for this information block is +``` +*ELEMENT, type=CPS4, ELSET=Surface1 +``` +The Abaqus element type `CPS4` corresponds to a quadrilateral element. +Each quadrilateral element in the unstructured mesh is dictated by four corner points with indexing +taken from the numbering given by the corner list above. +The elements connect a set of four corner points (starting from the bottom left) in an anti-clockwise fashion; +making the element *right-handed*. +This element handedness is indicated using the circular arrow in the figure above. +Just as with the corner list, the first integer in the element connectivity list indicates the element id number. +Thus, the element connectivity list for the three element example mesh is +``` +*ELEMENT, type=CPS4, ELSET=Surface1 +1, 5, 1, 4, 3 +2, 4, 2, 6, 3 +3, 7, 2, 4, 1 +``` + +#### Element neighbor connectivity + +The construction of the element neighbor ids and identifying physical boundary surfaces is done using functionality +directly from the [`p4est`](https://github.com/cburstedde/p4est) library. +For example, the neighbor connectivity is created in the mesh constructor using the wrapper `read_inp_p4est` function. + +#### HOHQMesh boundary information + +If present, any additional information in the mesh file that was created by `HOHQMesh` is prefaced with +`** ` to make it an Abaqus comment. +This ensures that the read in of the file by a standard Abaqus file parser, +as done in the `read_inp_p4est` function, is done correctly. + +The high-order, curved boundary information and labels of the physical boundary created by `HOHQMesh` +is found below the comment line +``` +** ***** HOHQMesh boundary information ***** ** +``` +Next comes the *polynomial degree* that the mesh will use to represent any curved sides +``` +** mesh polynomial degree = 8 +``` + +The mesh file then, again, provides the element connectivity as well as information for +curved surfaces either interior to the domain or along the physical boundaries. +A set of check digits are included directly below the four corner indexes to indicate whether +the local surface index (`-y`, `+x`, `+y`, or `-x`) within the element is +straight sided, `0`, or is curved, `1`. +If the local surface is straight sided no additional information is necessary during the mesh file read in. +But for any curved surfaces the mesh file provides `(x,y,z)` coordinate values in order to construct an +interpolant of this surface with the mesh polynomial order at the Chebyshev-Gauss-Lobatto +nodes. This list of `(x,y,z)` data will be given in the direction of the local coordinate system. +Given below is the element curvature information for the example mesh: +``` +** 5 1 4 3 +** 0 0 1 1 +** 1.000000000000000 1.000000000000000 0.0 +** 1.024948365654583 0.934461926834452 0.0 +** 1.116583018200151 0.777350964621867 0.0 +** 1.295753434047077 0.606254343587194 0.0 +** 1.537500000000000 0.462500000000000 0.0 +** 1.768263070247418 0.329729152118310 0.0 +** 1.920916981799849 0.185149035378133 0.0 +** 1.986035130050921 0.054554577460044 0.0 +** 2.000000000000000 0.0 0.0 +** 0.0 0.0 0.0 +** 0.035513826946206 0.105291711848750 0.0 +** 0.148591270347399 0.317731556850611 0.0 +** 0.340010713990041 0.452219430075470 0.0 +** 0.575000000000000 0.462500000000000 0.0 +** 0.788022294598950 0.483764065630034 0.0 +** 0.926408729652601 0.644768443149389 0.0 +** 0.986453164464803 0.883724792445746 0.0 +** 1.000000000000000 1.000000000000000 0.0 +** 4 2 6 3 +** 0 0 0 1 +** 2.000000000000000 0.0 0.0 +** 1.986035130050921 0.054554577460044 0.0 +** 1.920916981799849 0.185149035378133 0.0 +** 1.768263070247418 0.329729152118310 0.0 +** 1.537500000000000 0.462500000000000 0.0 +** 1.295753434047077 0.606254343587194 0.0 +** 1.116583018200151 0.777350964621867 0.0 +** 1.024948365654583 0.934461926834452 0.0 +** 1.000000000000000 1.000000000000000 0.0 +** 7 2 4 1 +** 0 0 0 0 +``` + +The last piece of information provided by `HOHQMesh` are labels for the different surfaces of an element. +These labels are useful to set boundary conditions along physical surfaces. +The labels can be short descriptive words up to 32 characters in length. +The label `---` indicates an internal surface where no boundary condition is required. + +It is important to note that these labels are given in the following order according to the +local surface index `-x` `+x` `-y` `+y` as required by the [`p4est`](https://github.com/cburstedde/p4est) library. +``` +** Bezier --- Slant --- +** --- Right --- Top +** Bottom --- Right --- +``` + +For completeness, we provide the entire Abaqus mesh file for the example mesh in the figure above: +``` +*Heading + File created by HOHQMesh +*NODE +1, 1.0, -1.0, 0.0 +2, 3.0, 0.0, 0.0 +3, 1.0, 1.0, 0.0 +4, 2.0, 0.0, 0.0 +5, 0.0, 0.0, 0.0 +6, 3.0, 1.0, 0.0 +7, 3.0, -1.0, 0.0 +*ELEMENT, type=CPS4, ELSET=Surface1 +1, 5, 1, 4, 3 +2, 4, 2, 6, 3 +3, 7, 2, 4, 1 +** ***** HOHQMesh boundary information ***** ** +** mesh polynomial degree = 8 +** 5 1 4 3 +** 0 0 1 1 +** 1.000000000000000 1.000000000000000 0.0 +** 1.024948365654583 0.934461926834452 0.0 +** 1.116583018200151 0.777350964621867 0.0 +** 1.295753434047077 0.606254343587194 0.0 +** 1.537500000000000 0.462500000000000 0.0 +** 1.768263070247418 0.329729152118310 0.0 +** 1.920916981799849 0.185149035378133 0.0 +** 1.986035130050921 0.054554577460044 0.0 +** 2.000000000000000 0.0 0.0 +** 0.0 0.0 0.0 +** 0.035513826946206 0.105291711848750 0.0 +** 0.148591270347399 0.317731556850611 0.0 +** 0.340010713990041 0.452219430075470 0.0 +** 0.575000000000000 0.462500000000000 0.0 +** 0.788022294598950 0.483764065630034 0.0 +** 0.926408729652601 0.644768443149389 0.0 +** 0.986453164464803 0.883724792445746 0.0 +** 1.000000000000000 1.000000000000000 0.0 +** 4 2 6 3 +** 0 0 0 1 +** 2.000000000000000 0.0 0.0 +** 1.986035130050921 0.054554577460044 0.0 +** 1.920916981799849 0.185149035378133 0.0 +** 1.768263070247418 0.329729152118310 0.0 +** 1.537500000000000 0.462500000000000 0.0 +** 1.295753434047077 0.606254343587194 0.0 +** 1.116583018200151 0.777350964621867 0.0 +** 1.024948365654583 0.934461926834452 0.0 +** 1.000000000000000 1.000000000000000 0.0 +** 7 2 4 1 +** 0 0 0 0 +** Bezier --- Slant --- +** --- Right --- Top +** Bottom --- Right --- +``` + +### Mesh in three spatial dimensions + +The 3D Abaqus file format with high-order boundary information from `HOHQMesh` is very similar to the +2D version discussed above. There are only three changes: + +1. The element connectivity would be given in terms of the eight corners that define a hexahedron. + The corners are numbered as shown in the figure below. The header of the element list changes to be + ``` + *ELEMENT, type=C3D8, ELSET=Volume1 + ``` + where `C3D8` corresponds to a Abaqus hexahedral element. +2. There are six check digits included directly below the eight corner indexes to indicate whether + the local face within the element is straight sided, `0`, or is curved, `1`. For curved faces + `(x,y,z)` coordinate values are available in order to construct an face interpolant with the mesh + polynomial order at the Chebyshev-Gauss-Lobatto nodes. +3. The boundary labels are given in the following order according to the local surface index + `-x` `+x` `-y` `+y` `-z` `+z` as required by the [`p4est`](https://github.com/cburstedde/p4est) library. + +For completeness, we also give a short description and derivation of the three-dimensional transfinite mapping +formulas used to compute the physical coordinates $\mathbf{x}=(x,y,z)$ of a (possibly curved) hexahedral element +give the reference coordinates $\boldsymbol{\xi} = (\xi, \eta, \zeta)$ which lie in $[-1,1]^3$. That is, we will +create an expression $\mathbf{x}= \mathbf{X}(\boldsymbol{\xi})$. + +Below we provide a sketch of a single hexahedral element with curved faces. This is done to introduce the numbering +conventions for corners, edges, and faces of the element. + +![abaqus-3dmesh-docs](https://user-images.githubusercontent.com/25242486/139839161-8c5f5979-2724-4cfb-9eac-6af58105ef12.png) + +When the hexahedron is a straight sided (linear) element we compute the transfinite mapping directly from the +element corner points according to +```math +\begin{aligned} +\mathbf{X}_{linear}(\boldsymbol{\xi}) &= \frac{1}{8}[\quad\, \mathbf{x}_1(1-\xi)(1-\eta)(1-\zeta) + + \mathbf{x}_2(1+\xi)(1-\eta)(1-\zeta)\\[-0.15cm] + & \qquad\; + \mathbf{x}_3(1+\xi)(1+\eta)(1-\zeta) + + \mathbf{x}_4(1-\xi)(1+\eta)(1-\zeta) \\ + & \qquad\; + \mathbf{x}_5(1-\xi)(1-\eta)(1+\zeta) + + \mathbf{x}_6(1+\xi)(1-\eta)(1+\zeta) \\ + & \qquad\; + \mathbf{x}_7(1+\xi)(1+\eta)(1+\zeta) + + \mathbf{x}_8(1-\xi)(1+\eta)(1+\zeta)\quad]. +\end{aligned} +``` + +Next, we create a transfinite mapping function, $\mathbf{X}(\boldsymbol{\xi})$, for a hexahedron that +has one or more curved faces. For this we assume that have a set of six interpolating polynomials +$\{\Gamma_i\}_{i=1}^6$ that approximate the faces. The interpolating polynomial for any curved faces is provided +by the information in a `HOHQMesh` Abaqus mesh file or is constructed on the fly via a +bi-linear interpolation routine for any linear faces. Explicitly, these six face interpolation polynomials depend +on the computational coordinates $\boldsymbol{\xi}$ as follows +```math + \begin{aligned} + \Gamma_1(\xi, \zeta), \quad && \quad \Gamma_3(\xi, \eta), \quad && \quad \Gamma_4(\eta, \zeta),\\[0.1cm] + \Gamma_2(\xi, \zeta), \quad && \quad \Gamma_5(\xi, \eta), \quad && \quad \Gamma_6(\eta, \zeta). + \end{aligned} +``` + +To determine the form of the mapping we first create linear interpolations between two opposing faces, e.g., $\Gamma_3$ and $\Gamma_5$ and sum them together to have +```math +\begin{aligned} + \boldsymbol\Sigma(\boldsymbol{\xi}) &= \frac{1}{2}[\quad\,(1-\xi)\Gamma_6(\eta,\zeta) + (1+\xi)\Gamma_4(\eta,\zeta) \\[-0.15cm] + &\qquad\;+ (1-\eta)\Gamma_1(\xi,\zeta) + (1+\eta)\Gamma_2(\xi,\zeta) \\%[-0.15cm] + &\qquad\; +(1-\zeta)\Gamma_3(\xi,\eta) + (1+\zeta)\Gamma_5(\xi,\eta)\quad]. +\end{aligned} +``` + +Unfortunately, the linear interpolations $\boldsymbol\Sigma(\boldsymbol{\xi})$ no longer match at the faces, e.g., evaluating at $\eta = -1$ we have +```math +\boldsymbol\Sigma(\xi,-1,\zeta) = \Gamma_1(\xi,\zeta) + \frac{1}{2}[\;(1-\xi)\Gamma_6(-1,\zeta) + (1+\xi)\Gamma_4(-1,\zeta) + +(1-\zeta)\Gamma_3(\xi,-1) + (1+\zeta)\Gamma_5(\xi,-1)\;], +``` +which is the desired face $\Gamma_1(\xi,\zeta)$ plus four edge error terms. +Analogous edge error terms occur at the other faces evaluating $\boldsymbol\Sigma(\boldsymbol{\xi})$ +at $\eta=1$, $\xi=\pm 1$, and $\zeta=\pm 1$. +In order to match the faces, we subtract a linear interpolant in the $\xi$, $\eta$, and $\zeta$ directions of the +edge error terms, e.g., the terms in braces in the above equation. So, continuing the example above, the correction term to be subtracted for face $\Gamma_1$ to match would be +```math +\left(\frac{1-\eta}{2}\right) \bigg[ \frac{1}{2} [ \; (1-\xi)\Gamma_6(-1,\zeta) + (1+\xi)\Gamma_4(-1,\zeta)+(1-\zeta)\Gamma_3(\xi,-1) + + (1+\zeta)\Gamma_5(\xi,-1)\;] \bigg]. +``` +For clarity, and to allow an easier comparison to the implementation, we introduce auxiliary notation for the 12 edge +values present in the complete correction term. That is, for given values of $\xi$, $\eta$, and $\zeta$ we have +```math + \begin{aligned} + \texttt{edge}_{1} &= \Gamma_1(\xi, -1), \quad && \quad \texttt{edge}_{5} = \Gamma_2(\xi, -1), \quad & \quad \texttt{edge}_{9} &= \Gamma_6(\eta, -1),\\[0.1cm] + \texttt{edge}_{2} &= \Gamma_1(1, \zeta), \quad && \quad\texttt{edge}_{6} = \Gamma_2(1, \zeta), \quad & \quad \texttt{edge}_{10} &= \Gamma_4(\eta, -1),\\[0.1cm] + \texttt{edge}_{3} &= \Gamma_1(\xi, 1), \quad && \quad \texttt{edge}_{7} = \Gamma_2(\xi, 1), \quad & \quad \texttt{edge}_{11} &= \Gamma_4(\eta, 1),\\[0.1cm] + \texttt{edge}_{4} &= \Gamma_1(-1, \zeta), \quad && \quad \texttt{edge}_{8} = \Gamma_2(-1, \zeta), \quad & \quad \texttt{edge}_{12} &= \Gamma_6(\eta, 1). + \end{aligned} +``` +With this notation for the edge terms (and after some algebraic manipulation) we write the complete edge correction term, +$\mathcal{C}_{\texttt{edge}}(\boldsymbol{\xi})$, as +```math +\begin{aligned} +\mathcal{C}_{\texttt{edge}}(\boldsymbol{\xi}) &= \frac{1}{4}[\quad\, (1-\eta)(1-\zeta)\texttt{edge}_{1}\\[-0.15cm] + & \qquad\; + (1+\xi)(1-\eta)\texttt{edge}_{2} \\ + & \qquad\; + (1-\eta)(1+\zeta)\texttt{edge}_{3} \\ + & \qquad\; + (1-\xi)(1-\eta)\texttt{edge}_{4} \\ + & \qquad\; + (1+\eta)(1-\zeta)\texttt{edge}_{5} \\ + & \qquad\; + (1+\xi)(1+\eta)\texttt{edge}_{6} \\ + & \qquad\; + (1+\eta)(1+\zeta)\texttt{edge}_{7} \\ + & \qquad\; + (1-\xi)(1+\eta)\texttt{edge}_{8} \\ + & \qquad\; + (1-\xi)(1-\zeta)\texttt{edge}_{9} \\ + & \qquad\; + (1+\xi)(1-\zeta)\texttt{edge}_{10} \\ + & \qquad\; + (1+\xi)(1+\zeta)\texttt{edge}_{11} \\ + & \qquad\; + (1-\xi)(1+\zeta)\texttt{edge}_{12}\quad]. +\end{aligned} +``` + +However, subtracting the edge correction terms $\mathcal{C}_{\texttt{edge}}(\boldsymbol{\xi})$ +from $\boldsymbol\Sigma(\boldsymbol{\xi})$ removes the interior element contributions twice. +Thus, to complete the construction of the transfinite mapping $\mathbf{X}(\boldsymbol{\xi})$ we add the +transfinite map of the straight sided hexahedral element to find +```math +\mathbf{X}(\boldsymbol{\xi}) = \boldsymbol\Sigma(\boldsymbol{\xi}) + - \mathcal{C}_{\texttt{edge}}(\boldsymbol{\xi}) + + \mathbf{X}_{linear}(\boldsymbol{\xi}). +``` \ No newline at end of file diff --git a/examples/p4est_2d_dgsem/elixir_euler_wall_bc_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_wall_bc_amr.jl new file mode 100644 index 00000000000..a95a4e18583 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_wall_bc_amr.jl @@ -0,0 +1,93 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler + +equations = CompressibleEulerEquations2D(1.4) + +@inline function uniform_flow_state(x, t, equations::CompressibleEulerEquations2D) + # set the freestream flow parameters + rho_freestream = 1.0 + u_freestream = 0.3 + p_freestream = inv(equations.gamma) + + theta = pi / 90.0 # analogous with a two degree angle of attack + si, co = sincos(theta) + v1 = u_freestream * co + v2 = u_freestream * si + + prim = SVector(rho_freestream, v1, v2, p_freestream) + return prim2cons(prim, equations) +end + +initial_condition = uniform_flow_state + +boundary_condition_uniform_flow = BoundaryConditionDirichlet(uniform_flow_state) +boundary_conditions = Dict( :Body => boundary_condition_uniform_flow, + :Button1 => boundary_condition_slip_wall, + :Button2 => boundary_condition_slip_wall, + :Eye1 => boundary_condition_slip_wall, + :Eye2 => boundary_condition_slip_wall, + :Smile => boundary_condition_slip_wall, + :Bowtie => boundary_condition_slip_wall ) + +volume_flux = flux_ranocha +solver = DGSEM(polydeg=5, surface_flux=flux_lax_friedrichs, + volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) + +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +default_mesh_file = joinpath(@__DIR__, "abaqus_gingerbread_man.inp") +isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/0e9e990a04b5105d1d2e3096a6e41272/raw/0d924b1d7e7d3cc1070a6cc22fe1d501687aa6dd/abaqus_gingerbread_man.inp", + default_mesh_file) +mesh_file = default_mesh_file + +mesh = P4estMesh{2}(mesh_file) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions=boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 3.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, + extra_analysis_integrals=(entropy,)) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=100, + save_initial_solution=true, + save_final_solution=true, + solution_variables=cons2prim) + +amr_indicator = IndicatorLöhner(semi, variable=density) + +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level=0, + med_level=1, med_threshold=0.05, + max_level=3, max_threshold=0.1) + +amr_callback = AMRCallback(semi, amr_controller, + interval=5, + adapt_initial_condition=true, + adapt_initial_condition_only_refine=true) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + amr_callback) + + +############################################################################### +# run the simulation +sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-7, reltol=1.0e-7, + save_everystep=false, callback=callbacks); +summary_callback() # print the timer summary diff --git a/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic_hohqmesh.jl b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic_hohqmesh.jl new file mode 100644 index 00000000000..ff101bea8aa --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic_hohqmesh.jl @@ -0,0 +1,69 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +initial_condition = initial_condition_convergence_test + +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = Dict( :Bottom => boundary_condition, + :Top => boundary_condition, + :Circle => boundary_condition, + :Cut => boundary_condition ) + +solver = DGSEM(polydeg=4, surface_flux=flux_lax_friedrichs) + +# Unstructured 3D half circle mesh from HOHQMesh +default_mesh_file = joinpath(@__DIR__, "abaqus_half_circle_3d.inp") +isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/11461efbfb02c42e06aca338b3d0b645/raw/81deeb1ebc4945952c30af5bb75fe222a18d975c/abaqus_half_circle_3d.inp", + default_mesh_file) +mesh_file = default_mesh_file + +mesh = P4estMesh{3}(mesh_file, initial_refinement_level=0) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms=source_terms_convergence_test, + boundary_conditions=boundary_conditions) + + +############################################################################### +# 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, + extra_analysis_integrals=(entropy,)) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=50, + save_initial_solution=true, + save_final_solution=true, + solution_variables=cons2prim) + +stepsize_callback = StepsizeCallback(cfl=1.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, # 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/src/auxiliary/p4est.jl b/src/auxiliary/p4est.jl index d3a47584f9a..afaf567972c 100644 --- a/src/auxiliary/p4est.jl +++ b/src/auxiliary/p4est.jl @@ -8,8 +8,8 @@ """ init_p4est() -Initialize p4est by calling `p4est_init` and setting the log level to `SC_LP_ERROR`. -This function will check if p4est is already initialized +Initialize `p4est` by calling `p4est_init` and setting the log level to `SC_LP_ERROR`. +This function will check if `p4est` is already initialized and if yes, do nothing, thus it is safe to call it multiple times. """ function init_p4est() @@ -17,7 +17,7 @@ function init_p4est() return nothing end - # Initialize p4est with log level ERROR to prevent a lot of output in AMR simulations + # Initialize `p4est` with log level ERROR to prevent a lot of output in AMR simulations p4est_init(C_NULL, SC_LP_ERROR) return nothing @@ -42,7 +42,7 @@ function unsafe_load_sc(::Type{T}, sc_array, i=1) where T end -# Create new p4est from a p4est_connectivity +# Create new `p4est` from a p4est_connectivity # 2D function new_p4est(conn::Ptr{p4est_connectivity_t}, initial_refinement_level) p4est_new_ext(0, # No MPI communicator @@ -61,7 +61,7 @@ function new_p4est(conn::Ptr{p8est_connectivity_t}, initial_refinement_level) end -# Save p4est data to file +# Save `p4est` data to file # 2D function save_p4est!(file, p4est::Ptr{p4est_t}) # Don't save user data of the quads @@ -75,7 +75,7 @@ function save_p4est!(file, p8est::Ptr{p8est_t}) end -# Load p4est from file +# Load `p4est` from file # 2D function load_p4est(file, ::Val{2}) conn_vec = Vector{Ptr{p4est_connectivity_t}}(undef, 1) @@ -89,28 +89,28 @@ function load_p4est(file, ::Val{3}) end -# Read p4est connectivity from Abaqus mesh file (.inp) +# Read `p4est` connectivity from Abaqus mesh file (.inp) # 2D read_inp_p4est(meshfile, ::Val{2}) = p4est_connectivity_read_inp(meshfile) # 3D read_inp_p4est(meshfile, ::Val{3}) = p8est_connectivity_read_inp(meshfile) -# Refine p4est if refine_fn_c returns 1 +# Refine `p4est` if refine_fn_c returns 1 # 2D refine_p4est!(p4est::Ptr{p4est_t}, recursive, refine_fn_c, init_fn_c) = p4est_refine(p4est, recursive, refine_fn_c, init_fn_c) # 3D refine_p4est!(p8est::Ptr{p8est_t}, recursive, refine_fn_c, init_fn_c) = p8est_refine(p8est, recursive, refine_fn_c, init_fn_c) -# Refine p4est if coarsen_fn_c returns 1 +# Refine `p4est` if coarsen_fn_c returns 1 # 2D coarsen_p4est!(p4est::Ptr{p4est_t}, recursive, coarsen_fn_c, init_fn_c) = p4est_coarsen(p4est, recursive, coarsen_fn_c, init_fn_c) # 3D coarsen_p4est!(p8est::Ptr{p8est_t}, recursive, coarsen_fn_c, init_fn_c) = p8est_coarsen(p8est, recursive, coarsen_fn_c, init_fn_c) -# Let p4est iterate over each cell volume and cell face. +# Let `p4est` iterate over each cell volume and cell face. # Call iter_volume_c for each cell and iter_face_c for each face. # 2D function iterate_p4est(p4est::Ptr{p4est_t}, user_data; diff --git a/src/meshes/face_interpolant.jl b/src/meshes/face_interpolant.jl new file mode 100644 index 00000000000..1d31de9152a --- /dev/null +++ b/src/meshes/face_interpolant.jl @@ -0,0 +1,52 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin + + +# CurvedFace{RealT<:Real} +# +# Contains the data needed to represent a curved face with data points (x,y,z) as a Lagrange polynomial +# interpolant written in barycentric form at a given set of nodes. +struct CurvedFace{RealT<:Real} + nodes ::Vector{RealT} + barycentric_weights ::Vector{RealT} + coordinates ::Array{RealT, 3} #[ndims, nnodes, nnodes] +end + + +# evalute the Gamma face interpolant at a particular point s = (s_1, s_2) and return the (x,y,z) coordinate +function evaluate_at(s, boundary_face::CurvedFace) + + @unpack nodes, barycentric_weights, coordinates = boundary_face + + x_coordinate_at_s_on_boundary_face = lagrange_interpolation_2d(s, nodes, view(coordinates, 1, :, :), + barycentric_weights) + y_coordinate_at_s_on_boundary_face = lagrange_interpolation_2d(s, nodes, view(coordinates, 2, :, :), + barycentric_weights) + z_coordinate_at_s_on_boundary_face = lagrange_interpolation_2d(s, nodes, view(coordinates, 3, :, :), + barycentric_weights) + + return x_coordinate_at_s_on_boundary_face, + y_coordinate_at_s_on_boundary_face, + z_coordinate_at_s_on_boundary_face +end + + +# Calculate a 2D Lagrange interpolating polynomial in barycentric 2 form +# of a function f(x,y) at a given coordinate (x,y) for a given node distribution. +function lagrange_interpolation_2d(x, nodes, function_values, barycentric_weights) + + f_intermediate = zeros(eltype(function_values), length(nodes)) + for j in eachindex(nodes) + f_intermediate[j] = lagrange_interpolation(x[2], nodes, view(function_values, j, :), + barycentric_weights) + end + point_value = lagrange_interpolation(x[1], nodes, f_intermediate, barycentric_weights) + + return point_value +end + + +end # @muladd diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index 3e30a1b514d..94d6d578c96 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -157,7 +157,7 @@ function save_mesh_file(mesh::P4estMesh, output_directory, timestep=0) p4est_file = joinpath(output_directory, p4est_filename) - # Save the complete connectivity/p4est data to disk. + # Save the complete connectivity and `p4est` data to disk. save_p4est!(p4est_file, mesh.p4est) # Open file (clobber existing content) @@ -264,7 +264,7 @@ function load_mesh_serial(mesh_file::AbstractString; n_cells_max, RealT) boundary_names = boundary_names_ .|> Symbol p4est_file = joinpath(dirname(mesh_file), p4est_filename) - # Prevent Julia crashes when p4est can't find the file + # Prevent Julia crashes when `p4est` can't find the file @assert isfile(p4est_file) p4est = load_p4est(p4est_file, Val(ndims)) diff --git a/src/meshes/meshes.jl b/src/meshes/meshes.jl index 25b8c74391b..a6dcbe132d8 100644 --- a/src/meshes/meshes.jl +++ b/src/meshes/meshes.jl @@ -9,6 +9,8 @@ include("tree_mesh.jl") include("structured_mesh.jl") include("surface_interpolant.jl") include("unstructured_mesh.jl") +include("face_interpolant.jl") +include("transfinite_mappings_3d.jl") include("p4est_mesh.jl") include("mesh_io.jl") include("dgmulti_meshes.jl") diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index b4bbb9e3159..f87302832c7 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -8,7 +8,7 @@ """ P4estMesh{NDIMS} <: AbstractMesh{NDIMS} -An unstructured curved mesh based on trees that uses the C library p4est +An unstructured curved mesh based on trees that uses the C library `p4est` to manage trees and mesh refinement. """ mutable struct P4estMesh{NDIMS, RealT<:Real, P, NDIMSP2, NNODES} <: AbstractMesh{NDIMS} @@ -32,7 +32,7 @@ mutable struct P4estMesh{NDIMS, RealT<:Real, P, NDIMSP2, NNODES} <: AbstractMesh mesh = new{NDIMS, eltype(tree_node_coordinates), typeof(p4est), NDIMS+2, length(nodes)}( p4est, tree_node_coordinates, nodes, boundary_names, current_filename, unsaved_changes) - # Destroy p4est structs when the mesh is garbage collected + # Destroy `p4est` structs when the mesh is garbage collected finalizer(destroy_mesh, mesh) return mesh @@ -247,15 +247,47 @@ end mapping=nothing, polydeg=1, RealT=Float64, initial_refinement_level=0, unsaved_changes=true) -Import an uncurved, unstructured, conforming mesh from an Abaqus mesh file (`.inp`), -map the mesh with the specified mapping, and create a `P4estMesh` from the curved mesh. - -Cells in the mesh file will be imported as trees in the `P4estMesh`. -The mesh will only have one boundary `:all`, as distinguishing different boundaries -is non-trivial. +Main mesh constructor for the `P4estMesh` that imports an unstructured, conforming +mesh from an Abaqus mesh file (`.inp`). Each element of the conforming mesh parsed +from the `meshfile` is created as a [`p4est`](https://github.com/cburstedde/p4est) +tree datatype. + +To create a curved unstructured mesh `P4estMesh` two strategies are available: + +- `p4est_mesh_from_hohqmesh_abaqus`: High-order, curved boundary information created by + [`HOHQMesh.jl`](https://github.com/trixi-framework/HOHQMesh.jl) is + available in the `meshfile`. The mesh polynomial degree `polydeg` + of the boundaries is provided from the `meshfile`. The computation of + the mapped tree coordinates is done with transfinite interpolation + with linear blending similar to [`UnstructuredMesh2D`](@ref). Boundary name + information is also parsed from the `meshfile` such that different boundary + conditions can be set at each named boundary on a given tree. +- `p4est_mesh_from_standard_abaqus`: By default, with `mapping=nothing` and `polydeg=1`, this creates a + straight-sided from the information parsed from the `meshfile`. If a mapping + function is specified then it computes the mapped tree coordinates via polynomial + interpolants with degree `polydeg`. The mesh created by this function will only + have one boundary `:all`, as distinguishing different physical boundaries is + non-trivial. + +Note that the `mapping` and `polydeg` keyword arguments are only used by the `p4est_mesh_from_standard_abaqus` +function. The `p4est_mesh_from_hohqmesh_abaqus` function obtains the mesh `polydeg` directly from the `meshfile` +and constructs the transfinite mapping internally. + +The particular strategy is selected according to the header present in the `meshfile` where +the constructor checks whether or not the `meshfile` was created with +[HOHQMesh.jl](https://github.com/trixi-framework/HOHQMesh.jl). +If the Abaqus file header is not present in the `meshfile` then the `P4estMesh` is created +with the function `p4est_mesh_from_standard_abaqus`. + +The default keyword argument `initial_refinement_level=0` corresponds to a forest +where the number of trees is the same as the number of elements in the original `meshfile`. +Increasing the `initial_refinement_level` allows one to uniformly refine the base mesh given +in the `meshfile` to create a forest with more trees before the simulation begins. +For example, if a two-dimensional base mesh contains 25 elements then setting +`initial_refinement_level=1` creates an initial forest of `2^2 * 25 = 100` trees. # Arguments -- `meshfile::String`: an uncurved Abaqus mesh file that can be imported by p4est. +- `meshfile::String`: an uncurved Abaqus mesh file that can be imported by `p4est`. - `mapping`: a function of `NDIMS` variables to describe the mapping that transforms the imported mesh to the physical domain. Use `nothing` for the identity map. - `polydeg::Integer`: polynomial degree used to store the geometry of the mesh. @@ -268,35 +300,118 @@ is non-trivial. - `unsaved_changes::Bool`: if set to `true`, the mesh will be saved to a mesh file. """ function P4estMesh{NDIMS}(meshfile::String; - mapping=nothing, polydeg=1, RealT=Float64, - initial_refinement_level=0, unsaved_changes=true) where NDIMS - # Prevent p4est from crashing Julia if the file doesn't exist + mapping=nothing, polydeg=1, RealT=Float64, + initial_refinement_level=0, unsaved_changes=true) where NDIMS + # Prevent `p4est` from crashing Julia if the file doesn't exist @assert isfile(meshfile) - conn = read_inp_p4est(meshfile, Val(NDIMS)) + # Read in the Header of the meshfile to determine which constructor is approproate + header = open(meshfile, "r") do io + readline(io) # *Header of the Abaqus file; discarded + readline(io) # Readin the actual header information + end + + # Check if the meshfile was generated using HOHQMesh + if header == " File created by HOHQMesh" + # Mesh curvature and boundary naming is handled with additional information available in meshfile + p4est, tree_node_coordinates, nodes, boundary_names = p4est_mesh_from_hohqmesh_abaqus(meshfile, initial_refinement_level, + NDIMS, RealT) + else + # Mesh curvature is handled directly by applying the mapping keyword argument + p4est, tree_node_coordinates, nodes, boundary_names = p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, + initial_refinement_level, + NDIMS, RealT) + end + + return P4estMesh{NDIMS}(p4est, tree_node_coordinates, nodes, + boundary_names, "", unsaved_changes) +end + + +# Create the mesh connectivity, mapped node coordinates within each tree, reference nodes in [-1,1] +# and a list of boundary names for the `P4estMesh`. High-order boundary curve information as well as +# the boundary names on each tree are provided by the `meshfile` created by +# [`HOHQMesh.jl`](https://github.com/trixi-framework/HOHQMesh.jl). +function p4est_mesh_from_hohqmesh_abaqus(meshfile, initial_refinement_level, n_dimensions, RealT) + # Create the mesh connectivity using `p4est` + conn = read_inp_p4est(meshfile, Val(n_dimensions)) # These need to be of the type Int for unsafe_wrap below to work n_trees::Int = conn.num_trees n_vertices::Int = conn.num_vertices - vertices = unsafe_wrap(Array, conn.vertices, (3, n_vertices)) - tree_to_vertex = unsafe_wrap(Array, conn.tree_to_vertex, (2^NDIMS, n_trees)) + # Extract a copy of the element vertices to compute the tree node coordinates + vertices = unsafe_wrap(Array, conn.vertices, (3, n_vertices)) + + # Readin all the information from the mesh file into a string array + file_lines = readlines(open(meshfile)) + + # Get the file index where the mesh polynomial degree is given in the meshfile + file_idx = findfirst(contains("** mesh polynomial degree"), file_lines) + + # Get the polynomial order of the mesh boundary information + current_line = split(file_lines[file_idx]) + mesh_polydeg = parse(Int, current_line[6]) + mesh_nnodes = mesh_polydeg + 1 + + # Create the Chebyshev-Gauss-Lobatto nodes used by HOHQMesh to represent the boundaries + cheby_nodes, _ = chebyshev_gauss_lobatto_nodes_weights(mesh_nnodes) + nodes = SVector{mesh_nnodes}(cheby_nodes) + + # Allocate the memory for the tree node coordinates + tree_node_coordinates = Array{RealT, n_dimensions+2}(undef, n_dimensions, + ntuple(_ -> length(nodes), n_dimensions)..., + n_trees) + + # Compute the tree node coordinates and return the updated file index + file_idx = calc_tree_node_coordinates!(tree_node_coordinates, file_lines, nodes, vertices, RealT) + + # Allocate the memory for the boundary labels + boundary_names = Array{Symbol}(undef, (2 * n_dimensions, n_trees)) + + # Read in the boundary names from the last portion of the meshfile + # Note here the boundary names where "---" means an internal connection + for tree in 1:n_trees + current_line = split(file_lines[file_idx]) + boundary_names[:, tree] = map(Symbol, current_line[2:end]) + file_idx += 1 + end + + p4est = new_p4est(conn, initial_refinement_level) + + return p4est, tree_node_coordinates, nodes, boundary_names +end + + +# Create the mesh connectivity, mapped node coordinates within each tree, reference nodes in [-1,1] +# and a list of boundary names for the `P4estMesh`. The tree node coordinates are computed according to +# the `mapping` passed to this function using polynomial interpolants of degree `polydeg`. All boundary +# names are given the name `:all`. +function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, initial_refinement_level, n_dimensions, RealT) + # Create the mesh connectivity using `p4est` + conn = read_inp_p4est(meshfile, Val(n_dimensions)) + + # These need to be of the type Int for unsafe_wrap below to work + n_trees::Int = conn.num_trees + n_vertices::Int = conn.num_vertices + + vertices = unsafe_wrap(Array, conn.vertices, (3, n_vertices)) + tree_to_vertex = unsafe_wrap(Array, conn.tree_to_vertex, (2^n_dimensions, n_trees)) basis = LobattoLegendreBasis(RealT, polydeg) nodes = basis.nodes - tree_node_coordinates = Array{RealT, NDIMS+2}(undef, NDIMS, - ntuple(_ -> length(nodes), NDIMS)..., - n_trees) + tree_node_coordinates = Array{RealT, n_dimensions+2}(undef, n_dimensions, + ntuple(_ -> length(nodes), n_dimensions)..., + n_trees) calc_tree_node_coordinates!(tree_node_coordinates, nodes, mapping, vertices, tree_to_vertex) p4est = new_p4est(conn, initial_refinement_level) # There's no simple and generic way to distinguish boundaries. Name all of them :all. - boundary_names = fill(:all, 2 * NDIMS, n_trees) + boundary_names = fill(:all, 2 * n_dimensions, n_trees) - return P4estMesh{NDIMS}(p4est, tree_node_coordinates, nodes, - boundary_names, "", unsaved_changes) + return p4est, tree_node_coordinates, nodes, boundary_names end @@ -359,9 +474,9 @@ end function connectivity_structured(n_cells_x, n_cells_y, periodicity) linear_indices = LinearIndices((n_cells_x, n_cells_y)) - # Vertices represent the coordinates of the forest. This is used by p4est + # Vertices represent the coordinates of the forest. This is used by `p4est` # to write VTK files. - # Trixi doesn't use p4est's coordinates, so the vertices can be empty. + # Trixi doesn't use the coordinates from `p4est`, so the vertices can be empty. n_vertices = 0 n_trees = n_cells_x * n_cells_y # No corner connectivity is needed @@ -375,7 +490,7 @@ function connectivity_structured(n_cells_x, n_cells_y, periodicity) for cell_y in 1:n_cells_y, cell_x in 1:n_cells_x tree = linear_indices[cell_x, cell_y] - # Subtract 1 because p4est uses zero-based indexing + # Subtract 1 because `p4est` uses zero-based indexing # Negative x-direction if cell_x > 1 tree_to_tree[1, tree] = linear_indices[cell_x - 1, cell_y] - 1 @@ -426,7 +541,7 @@ function connectivity_structured(n_cells_x, n_cells_y, periodicity) end tree_to_corner = C_NULL - # p4est docs: "in trivial cases it is just a pointer to a p4est_topix value of 0." + # `p4est` docs: "in trivial cases it is just a pointer to a p4est_topix value of 0." # We don't need corner connectivity, so this is a trivial case. ctt_offset = zeros(p4est_topidx_t, 1) @@ -448,9 +563,9 @@ end function connectivity_structured(n_cells_x, n_cells_y, n_cells_z, periodicity) linear_indices = LinearIndices((n_cells_x, n_cells_y, n_cells_z)) - # Vertices represent the coordinates of the forest. This is used by p4est + # Vertices represent the coordinates of the forest. This is used by `p4est` # to write VTK files. - # Trixi doesn't use p4est's coordinates, so the vertices can be empty. + # Trixi doesn't use the coordinates from `p4est`, so the vertices can be empty. n_vertices = 0 n_trees = n_cells_x * n_cells_y * n_cells_z # No edge connectivity is needed @@ -466,7 +581,7 @@ function connectivity_structured(n_cells_x, n_cells_y, n_cells_z, periodicity) for cell_z in 1:n_cells_z, cell_y in 1:n_cells_y, cell_x in 1:n_cells_x tree = linear_indices[cell_x, cell_y, cell_z] - # Subtract 1 because p4est uses zero-based indexing + # Subtract 1 because `p4est` uses zero-based indexing # Negative x-direction if cell_x > 1 tree_to_tree[1, tree] = linear_indices[cell_x - 1, cell_y, cell_z] - 1 @@ -541,14 +656,14 @@ function connectivity_structured(n_cells_x, n_cells_y, n_cells_z, periodicity) end tree_to_edge = C_NULL - # p4est docs: "in trivial cases it is just a pointer to a p4est_topix value of 0." + # `p4est` docs: "in trivial cases it is just a pointer to a p4est_topix value of 0." # We don't need edge connectivity, so this is a trivial case. ett_offset = zeros(p4est_topidx_t, 1) edge_to_tree = C_NULL edge_to_edge = C_NULL tree_to_corner = C_NULL - # p4est docs: "in trivial cases it is just a pointer to a p4est_topix value of 0." + # `p4est` docs: "in trivial cases it is just a pointer to a p4est_topix value of 0." # We don't need corner connectivity, so this is a trivial case. ctt_offset = zeros(p4est_topidx_t, 1) @@ -575,9 +690,9 @@ function connectivity_cubed_sphere(trees_per_face_dimension, layers) linear_indices = LinearIndices((trees_per_face_dimension, trees_per_face_dimension, layers, 6)) - # Vertices represent the coordinates of the forest. This is used by p4est + # Vertices represent the coordinates of the forest. This is used by `p4est` # to write VTK files. - # Trixi doesn't use p4est's coordinates, so the vertices can be empty. + # Trixi doesn't use the coordinates from `p4est`, so the vertices can be empty. n_vertices = 0 n_trees = 6 * n_cells_x * n_cells_y * n_cells_z # No edge connectivity is needed @@ -635,7 +750,7 @@ function connectivity_cubed_sphere(trees_per_face_dimension, layers) for cell_z in 1:n_cells_z, cell_y in 1:n_cells_y, cell_x in 1:n_cells_x tree = linear_indices[cell_x, cell_y, cell_z, direction] - # Subtract 1 because p4est uses zero-based indexing + # Subtract 1 because `p4est` uses zero-based indexing # Negative x-direction if cell_x > 1 # Connect to tree at the same face tree_to_tree[1, tree] = linear_indices[cell_x - 1, cell_y, cell_z, direction] - 1 @@ -777,14 +892,14 @@ function connectivity_cubed_sphere(trees_per_face_dimension, layers) end tree_to_edge = C_NULL - # p4est docs: "in trivial cases it is just a pointer to a p4est_topix value of 0." + # `p4est` docs: "in trivial cases it is just a pointer to a p4est_topix value of 0." # We don't need edge connectivity, so this is a trivial case. ett_offset = zeros(p4est_topidx_t, 1) edge_to_tree = C_NULL edge_to_edge = C_NULL tree_to_corner = C_NULL - # p4est docs: "in trivial cases it is just a pointer to a p4est_topix value of 0." + # `p4est` docs: "in trivial cases it is just a pointer to a p4est_topix value of 0." # We don't need corner connectivity, so this is a trivial case. ctt_offset = zeros(p4est_topidx_t, 1) @@ -1021,9 +1136,262 @@ function cubed_sphere_mapping(xi, eta, zeta, inner_radius, thickness, direction) end +# Calculate physical coordinates of each element of an unstructured mesh read +# in from a HOHQMesh file. This calculation is done with the transfinite interpolation +# routines found in `mappings_geometry_curved_2d.jl` or `mappings_geometry_straight_2d.jl` +function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, + file_lines::Vector{String}, nodes, vertices, RealT) + # Get the number of trees and the number of interpolation nodes + n_trees = last(size(node_coordinates)) + nnodes = length(nodes) + + # Setup the starting file index to read in element indices and the additional + # curved boundary information provided by HOHQMesh. + file_idx = findfirst(contains("** mesh polynomial degree"), file_lines) + 1 + + # Create a work set of Gamma curves to create the node coordinates + CurvedSurfaceT = CurvedSurface{RealT} + surface_curves = Array{CurvedSurfaceT}(undef, 4) + + # Create other work arrays to perform the mesh construction + element_node_ids = Array{Int}(undef, 4) + curved_check = Vector{Int}(undef, 4) + quad_vertices = Array{RealT}(undef, (4, 2)) + quad_vertices_flipped = Array{RealT}(undef, (4, 2)) + curve_values = Array{RealT}(undef, (nnodes, 2)) + + # Create the barycentric weights used for the surface interpolations + bary_weights_ = barycentric_weights(nodes) + bary_weights = SVector{nnodes}(bary_weights_) + + # Loop through all the trees, i.e., the elements generated by HOHQMesh and create the node coordinates. + # When we extract information from the `current_line` we start at index 2 in order to + # avoid the Abaqus comment character "** " + for tree in 1:n_trees + # Pull the vertex node IDs + current_line = split(file_lines[file_idx]) + element_node_ids[1] = parse(Int, current_line[2]) + element_node_ids[2] = parse(Int, current_line[3]) + element_node_ids[3] = parse(Int, current_line[4]) + element_node_ids[4] = parse(Int, current_line[5]) + + # Pull the (x,y) values of the four vertices of the current tree out of the global vertices array + for i in 1:4 + quad_vertices[i, :] .= vertices[1:2, element_node_ids[i]] + end + # Pull the information to check if boundary is curved in order to read in additional data + file_idx += 1 + current_line = split(file_lines[file_idx]) + curved_check[1] = parse(Int, current_line[2]) + curved_check[2] = parse(Int, current_line[3]) + curved_check[3] = parse(Int, current_line[4]) + curved_check[4] = parse(Int, current_line[5]) + if sum(curved_check) == 0 + # Create the node coordinates on this particular element + calc_node_coordinates!(node_coordinates, tree, nodes, quad_vertices) + else + # Quadrilateral element has at least one curved side + # Flip node ordering to make sure the element is right-handed for the interpolations + m1 = 1 + m2 = 2 + @views quad_vertices_flipped[1, :] .= quad_vertices[4, :] + @views quad_vertices_flipped[2, :] .= quad_vertices[2, :] + @views quad_vertices_flipped[3, :] .= quad_vertices[3, :] + @views quad_vertices_flipped[4, :] .= quad_vertices[1, :] + for i in 1:4 + if curved_check[i] == 0 + # When curved_check[i] is 0 then the "curve" from vertex `i` to vertex `i+1` is a straight line. + # Evaluate a linear interpolant between the two points at each of the nodes. + for k in 1:nnodes + curve_values[k, 1] = linear_interpolate(nodes[k], quad_vertices_flipped[m1, 1], quad_vertices_flipped[m2, 1]) + curve_values[k, 2] = linear_interpolate(nodes[k], quad_vertices_flipped[m1, 2], quad_vertices_flipped[m2, 2]) + end + else + # When curved_check[i] is 1 this curved boundary information is supplied by the mesh + # generator. So we just read it into a work array + for k in 1:nnodes + file_idx += 1 + current_line = split(file_lines[file_idx]) + curve_values[k, 1] = parse(RealT,current_line[2]) + curve_values[k, 2] = parse(RealT,current_line[3]) + end + end + # Construct the curve interpolant for the current side + surface_curves[i] = CurvedSurfaceT(nodes, bary_weights, copy(curve_values)) + # Indexing update that contains a "flip" to ensure correct element orientation. + # If we need to construct the straight line "curves" when curved_check[i] == 0 + m1 += 1 + if i == 3 + m2 = 1 + else + m2 += 1 + end + end + # Create the node coordinates on this particular element + calc_node_coordinates!(node_coordinates, tree, nodes, surface_curves) + end + # Move file index to the next tree + file_idx += 1 + end + + return file_idx +end + + +# Calculate physical coordinates of each element of an unstructured mesh read +# in from a HOHQMesh file. This calculation is done with the transfinite interpolation +# routines found in `transfinite_mappings_3d.jl` +function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, + file_lines::Vector{String}, nodes, vertices, RealT) + # Get the number of trees and the number of interpolation nodes + n_trees = last(size(node_coordinates)) + nnodes = length(nodes) + + # Setup the starting file index to read in element indices and the additional + # curved boundary information provided by HOHQMesh. + file_idx = findfirst(contains("** mesh polynomial degree"), file_lines) + 1 + + # Create a work set of Gamma curves to create the node coordinates + CurvedFaceT = CurvedFace{RealT} + face_curves = Array{CurvedFaceT}(undef, 6) + + # Create other work arrays to perform the mesh construction + element_node_ids = Array{Int}(undef, 8) + curved_check = Vector{Int}(undef, 6) + hex_vertices = Array{RealT}(undef, (3, 8)) + face_vertices = Array{RealT}(undef, (3, 4)) + curve_values = Array{RealT}(undef, (3, nnodes, nnodes)) + + # Create the barycentric weights used for the surface interpolations + bary_weights_ = barycentric_weights(nodes) + bary_weights = SVector{nnodes}(bary_weights_) + + # Loop through all the trees, i.e., the elements generated by HOHQMesh and create the node coordinates. + # When we extract information from the `current_line` we start at index 2 in order to + # avoid the Abaqus comment character "** " + for tree in 1:n_trees + # pull the vertex node IDs + current_line = split(file_lines[file_idx]) + element_node_ids[1] = parse(Int, current_line[2]) + element_node_ids[2] = parse(Int, current_line[3]) + element_node_ids[3] = parse(Int, current_line[4]) + element_node_ids[4] = parse(Int, current_line[5]) + element_node_ids[5] = parse(Int, current_line[6]) + element_node_ids[6] = parse(Int, current_line[7]) + element_node_ids[7] = parse(Int, current_line[8]) + element_node_ids[8] = parse(Int, current_line[9]) + + # Pull the (x, y, z) values of the eight vertices of the current tree out of the global vertices array + for i in 1:8 + hex_vertices[:, i] .= vertices[:, element_node_ids[i]] + end + # Pull the information to check if boundary is curved in order to read in additional data + file_idx += 1 + current_line = split(file_lines[file_idx]) + curved_check[1] = parse(Int, current_line[2]) + curved_check[2] = parse(Int, current_line[3]) + curved_check[3] = parse(Int, current_line[4]) + curved_check[4] = parse(Int, current_line[5]) + curved_check[5] = parse(Int, current_line[6]) + curved_check[6] = parse(Int, current_line[7]) + if sum(curved_check) == 0 + # Create the node coordinates on this element + calc_node_coordinates!(node_coordinates, tree, nodes, hex_vertices) + else + # Hexahedral element has at least one curved side + for face in 1:6 + if curved_check[face] == 0 + # Face is a flat plane. Evaluate a bilinear interpolant between the four vertices of the face at each of the nodes. + get_vertices_for_bilinear_interpolant!(face_vertices, face, hex_vertices) + for q in 1:nnodes, p in 1:nnodes + @views bilinear_interpolation!(curve_values[:, p, q], face_vertices, nodes[p], nodes[q]) + end + else # curved_check[face] == 1 + # Curved face boundary information is supplied by the mesh file. Just read it into a work array + for q in 1:nnodes, p in 1:nnodes + file_idx += 1 + current_line = split(file_lines[file_idx]) + curve_values[1, p, q] = parse(RealT,current_line[2]) + curve_values[2, p, q] = parse(RealT,current_line[3]) + curve_values[3, p, q] = parse(RealT,current_line[4]) + end + end + # Construct the curve interpolant for the current side + face_curves[face] = CurvedFaceT(nodes, bary_weights, copy(curve_values)) + end + # Create the node coordinates on this particular element + calc_node_coordinates!(node_coordinates, tree, nodes, face_curves) + end + # Move file index to the next tree + file_idx += 1 + end + + return file_idx +end + + +# Given the eight `hex_vertices` for a hexahedral element extract +# the four `face_vertices` for a particular `face_index`. +function get_vertices_for_bilinear_interpolant!(face_vertices, face_index, hex_vertices) + if face_index == 1 + @views face_vertices[:, 1] .= hex_vertices[:, 1] + @views face_vertices[:, 2] .= hex_vertices[:, 2] + @views face_vertices[:, 3] .= hex_vertices[:, 6] + @views face_vertices[:, 4] .= hex_vertices[:, 5] + elseif face_index == 2 + @views face_vertices[:, 1] .= hex_vertices[:, 4] + @views face_vertices[:, 2] .= hex_vertices[:, 3] + @views face_vertices[:, 3] .= hex_vertices[:, 7] + @views face_vertices[:, 4] .= hex_vertices[:, 8] + elseif face_index == 3 + @views face_vertices[:, 1] .= hex_vertices[:, 1] + @views face_vertices[:, 2] .= hex_vertices[:, 2] + @views face_vertices[:, 3] .= hex_vertices[:, 3] + @views face_vertices[:, 4] .= hex_vertices[:, 4] + elseif face_index == 4 + @views face_vertices[:, 1] .= hex_vertices[:, 2] + @views face_vertices[:, 2] .= hex_vertices[:, 3] + @views face_vertices[:, 3] .= hex_vertices[:, 6] + @views face_vertices[:, 4] .= hex_vertices[:, 7] + elseif face_index == 5 + @views face_vertices[:, 1] .= hex_vertices[:, 5] + @views face_vertices[:, 2] .= hex_vertices[:, 6] + @views face_vertices[:, 3] .= hex_vertices[:, 7] + @views face_vertices[:, 4] .= hex_vertices[:, 8] + else # face_index == 6 + @views face_vertices[:, 1] .= hex_vertices[:, 1] + @views face_vertices[:, 2] .= hex_vertices[:, 4] + @views face_vertices[:, 3] .= hex_vertices[:, 8] + @views face_vertices[:, 4] .= hex_vertices[:, 5] + end +end + + +# Evaluate a bilinear interpolant at a point (u,v) given the four vertices where the face is right-handed +# 4 3 +# o----------------o +# | | +# | | +# | | +# | | +# | | +# | | +# o----------------o +# 1 2 +# and return the 3D coordinate point (x, y, z) +function bilinear_interpolation!(coordinate, face_vertices, u, v) + for j in 1:3 + coordinate[j] = 0.25 * ( face_vertices[j,1] * (1 - u) * (1 - v) + + face_vertices[j,2] * (1 + u) * (1 - v) + + face_vertices[j,3] * (1 + u) * (1 + v) + + face_vertices[j,4] * (1 - u) * (1 + v) ) + end +end + + function balance!(mesh::P4estMesh{2}, init_fn=C_NULL) p4est_balance(mesh.p4est, P4EST_CONNECT_FACE, init_fn) - # Due to a bug in p4est, the forest needs to be rebalanced twice sometimes + # Due to a bug in `p4est`, the forest needs to be rebalanced twice sometimes # See https://github.com/cburstedde/p4est/issues/112 p4est_balance(mesh.p4est, P4EST_CONNECT_FACE, init_fn) end @@ -1095,7 +1463,7 @@ function coarsen_fn(p4est, which_tree, quadrants_ptr) # Load controller value from quadrant's user data ([global quad ID, controller_value]). controller_value(i) = unsafe_load(Ptr{Int}(quadrants[i].p.user_data), 2) - # p4est calls this function for each 2^ndims quads that could be coarsened to a single one. + # `p4est` calls this function for each 2^ndims quads that could be coarsened to a single one. # Only coarsen if all these 2^ndims quads have been marked for coarsening. if all(i -> controller_value(i) < 0, eachindex(quadrants)) # return true (coarsen) diff --git a/src/meshes/transfinite_mappings_3d.jl b/src/meshes/transfinite_mappings_3d.jl new file mode 100644 index 00000000000..b74a72be8e7 --- /dev/null +++ b/src/meshes/transfinite_mappings_3d.jl @@ -0,0 +1,172 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin + +# Illustration of the corner (circled), edge (braces), and face index numbering convention +# used in these functions. +# +# ⑧────────────────────────{7}────────────────────────⑦ +# ╱│ ╱│ +# ╱ │ ╱ │ +# ╱ │ ╱ │ +# ╱ │ 5 (+z) ╱ │ +# ╱ │ ╱ │ +# ╱ │ ╱ │ +# {12} │ {11} │ +# ╱ │ ╱ │ +# ╱ │ ╱ │ +# ╱ │ 2 (+y) ╱ │ +# ╱ │ ╱ │ +# ╱ {8} ╱ {6} +# ╱ │ ╱ │ +# ⑤─────────────────────────{3}───────────────────────⑥ 4 (+x) │ +# │ │ │ │ +# │ │ │ │ +# │ │ │ │ +# │ 6 (-x) │ │ │ +# │ │ │ │ +# │ │ │ │ +# │ │ │ │ +# │ │ │ │ +# │ ④────────────────────────{5}──────────│─────────────③ +# │ ╱ │ ╱ +# │ ╱ 1 (-y) │ ╱ +# {4} ╱ {2} ╱ +# │ ╱ │ ╱ +# │ ╱ │ ╱ +# │ {9} │ {10} +# │ ╱ │ ╱ +# │ ╱ │ ╱ Global coordinates: +# │ ╱ │ ╱ z +# │ ╱ 3 (-z) │ ╱ ↑ y +# │ ╱ │ ╱ │ ╱ +# │ ╱ │ ╱ │ ╱ +# │╱ │╱ └─────> x +# ①───────────────────────{1}─────────────────────────② + + +# Transfinite mapping formula from a point (xi, eta, zeta) in reference space [-1,1]^3 to a +# physical coordinate (x, y, z) for a hexahedral element with straight sides +function straight_side_hex_map(xi, eta, zeta, corner_points) + + coordinate = zeros(eltype(xi), 3) + for j in 1:3 + coordinate[j] += (0.125 * ( corner_points[j, 1] * (1 - xi) * (1 - eta) * (1 - zeta) + + corner_points[j, 2] * (1 + xi) * (1 - eta) * (1 - zeta) + + corner_points[j, 3] * (1 + xi) * (1 + eta) * (1 - zeta) + + corner_points[j, 4] * (1 - xi) * (1 + eta) * (1 - zeta) + + corner_points[j, 5] * (1 - xi) * (1 - eta) * (1 + zeta) + + corner_points[j, 6] * (1 + xi) * (1 - eta) * (1 + zeta) + + corner_points[j, 7] * (1 + xi) * (1 + eta) * (1 + zeta) + + corner_points[j, 8] * (1 - xi) * (1 + eta) * (1 + zeta) ) ) + end + + return coordinate +end + + +# Construct the (x, y, z) node coordinates in the volume of a straight sided hexahedral element +function calc_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, element, nodes, corners) + + for k in eachindex(nodes), j in eachindex(nodes), i in eachindex(nodes) + node_coordinates[:, i, j, k, element] .= straight_side_hex_map(nodes[i], nodes[j], nodes[k], corners) + end + + return node_coordinates +end + + +# Transfinite mapping formula from a point (xi, eta, zeta) in reference space [-1,1]^3 to a point +# (x,y,z) in physical coordinate space for a hexahedral element with general curved sides +# See Section 4.3 +# - Andrew R. Winters (2014) +# Discontinous Galerkin spectral element approximations for the reflection and +# transmission of waves from moving material interfaces +# [PhD thesis, Florida State University](https://diginole.lib.fsu.edu/islandora/object/fsu%3A185342) +function transfinite_hex_map(xi, eta, zeta, face_curves::AbstractVector{<:CurvedFace}) + + coordinate = zeros(eltype(xi), 3) + face_values = zeros(eltype(xi), (3, 6)) + edge_values = zeros(eltype(xi), (3, 12)) + corners = zeros(eltype(xi), (3, 8)) + + # Compute values along the face edges + edge_values[:, 1] .= evaluate_at(SVector(xi, -1), face_curves[1]) + edge_values[:, 2] .= evaluate_at(SVector( 1, zeta), face_curves[1]) + edge_values[:, 3] .= evaluate_at(SVector(xi, 1), face_curves[1]) + edge_values[:, 4] .= evaluate_at(SVector(-1, zeta), face_curves[1]) + + edge_values[:, 5] .= evaluate_at(SVector(xi, -1), face_curves[2]) + edge_values[:, 6] .= evaluate_at(SVector( 1, zeta), face_curves[2]) + edge_values[:, 7] .= evaluate_at(SVector(xi, 1), face_curves[2]) + edge_values[:, 8] .= evaluate_at(SVector(-1, zeta), face_curves[2]) + + edge_values[:, 9] .= evaluate_at(SVector(eta, -1), face_curves[6]) + edge_values[:, 10] .= evaluate_at(SVector(eta, -1), face_curves[4]) + edge_values[:, 11] .= evaluate_at(SVector(eta, 1), face_curves[4]) + edge_values[:, 12] .= evaluate_at(SVector(eta, 1), face_curves[6]) + + # Compute values on the face + face_values[:, 1] .= evaluate_at(SVector( xi, zeta), face_curves[1]) + face_values[:, 2] .= evaluate_at(SVector( xi, zeta), face_curves[2]) + face_values[:, 3] .= evaluate_at(SVector( xi, eta), face_curves[3]) + face_values[:, 4] .= evaluate_at(SVector(eta, zeta), face_curves[4]) + face_values[:, 5] .= evaluate_at(SVector( xi, eta), face_curves[5]) + face_values[:, 6] .= evaluate_at(SVector(eta, zeta), face_curves[6]) + + # Pull the eight corner values and compute the straight sided hex mapping + corners[:,1] .= face_curves[1].coordinates[:, 1, 1] + corners[:,2] .= face_curves[1].coordinates[:, end, 1] + corners[:,3] .= face_curves[2].coordinates[:, end, 1] + corners[:,4] .= face_curves[2].coordinates[:, 1, 1] + corners[:,5] .= face_curves[1].coordinates[:, 1, end] + corners[:,6] .= face_curves[1].coordinates[:, end, end] + corners[:,7] .= face_curves[2].coordinates[:, end, end] + corners[:,8] .= face_curves[2].coordinates[:, 1, end] + + coordinate_straight = straight_side_hex_map(xi, eta, zeta, corners) + + # Compute the transfinite mapping + for j in 1:3 + # Linear interpolation between opposite faces + coordinate[j] = ( 0.5 * ( face_values[j, 6] * (1 - xi ) + face_values[j, 4] * (1 + xi ) + + face_values[j, 1] * (1 - eta ) + face_values[j, 2] * (1 + eta ) + + face_values[j, 3] * (1 - zeta) + face_values[j, 5] * (1 + zeta) ) ) + + # Edge corrections to ensure faces match + coordinate[j] -= ( 0.25 * ( edge_values[j, 1 ] * (1 - eta) * (1 - zeta) + + edge_values[j, 2 ] * (1 + xi ) * (1 - eta ) + + edge_values[j, 3 ] * (1 - eta) * (1 + zeta) + + edge_values[j, 4 ] * (1 - xi ) * (1 - eta ) + + edge_values[j, 5 ] * (1 + eta) * (1 - zeta) + + edge_values[j, 6 ] * (1 + xi ) * (1 + eta ) + + edge_values[j, 7 ] * (1 + eta) * (1 + zeta) + + edge_values[j, 8 ] * (1 - xi ) * (1 + eta ) + + edge_values[j, 9 ] * (1 - xi ) * (1 - zeta) + + edge_values[j, 10] * (1 + xi ) * (1 - zeta) + + edge_values[j, 11] * (1 + xi ) * (1 + zeta) + + edge_values[j, 12] * (1 - xi ) * (1 + zeta) ) ) + + # Subtracted interior twice, so add back the straight-sided hexahedral mapping + coordinate[j] += coordinate_straight[j] + end + + return coordinate +end + + +# Construct the (x, y, z) node coordinates in the volume of a curved sided hexahedral element +function calc_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, element, nodes, + face_curves::AbstractVector{<:CurvedFace}) + + for k in eachindex(nodes), j in eachindex(nodes), i in eachindex(nodes) + node_coordinates[:, i, j, k, element] .= transfinite_hex_map(nodes[i], nodes[j], nodes[k], face_curves) + end + + return node_coordinates +end + + +end # @muladd diff --git a/src/meshes/unstructured_mesh.jl b/src/meshes/unstructured_mesh.jl index 718c9929f7a..202abe8079b 100644 --- a/src/meshes/unstructured_mesh.jl +++ b/src/meshes/unstructured_mesh.jl @@ -45,9 +45,9 @@ function UnstructuredMesh2D(filename; RealT=Float64, periodicity=false, unsaved_ # readin the number of nodes, number of interfaces, number of elements and local polynomial degree current_line = split(file_lines[2]) - n_corners = parse(Int, current_line[1]) - n_surfaces = parse(Int, current_line[2]) - n_elements = parse(Int, current_line[3]) + n_corners = parse(Int, current_line[1]) + n_surfaces = parse(Int, current_line[2]) + n_elements = parse(Int, current_line[3]) mesh_polydeg = parse(Int, current_line[4]) mesh_nnodes = mesh_polydeg + 1 @@ -56,17 +56,17 @@ function UnstructuredMesh2D(filename; RealT=Float64, periodicity=false, unsaved_ # the mesh file. Thus, this cannot be type stable at all. Hence, we allocate # the memory now and introduce a function barrier before continuing to read # data from the file. - corner_nodes = Array{RealT}(undef, (2, n_corners)) - interface_info = Array{Int}(undef, (6, n_surfaces)) - element_node_ids = Array{Int}(undef, (4, n_elements)) - curved_check = Vector{Int}(undef, 4) - cornerNodeVals = Array{RealT}(undef, (4, 2)) - tempNodes = Array{RealT}(undef, (4, 2)) - curve_vals = Array{RealT}(undef, (mesh_nnodes, 2)) + corner_nodes = Array{RealT}(undef, (2, n_corners)) + interface_info = Array{Int}(undef, (6, n_surfaces)) + element_node_ids = Array{Int}(undef, (4, n_elements)) + curved_check = Vector{Int}(undef, 4) + quad_corners = Array{RealT}(undef, (4, 2)) + quad_corners_flipped = Array{RealT}(undef, (4, 2)) + curve_values = Array{RealT}(undef, (mesh_nnodes, 2)) element_is_curved = Array{Bool}(undef, n_elements) - CurvedSurfaceT = CurvedSurface{RealT} - surface_curves = Array{CurvedSurfaceT}(undef, (4, n_elements)) - boundary_names = Array{Symbol}(undef, (4, n_elements)) + CurvedSurfaceT = CurvedSurface{RealT} + surface_curves = Array{CurvedSurfaceT}(undef, (4, n_elements)) + boundary_names = Array{Symbol}(undef, (4, n_elements)) # create the Chebyshev-Gauss-Lobatto nodes used to represent any curved boundaries that are # required to construct the sides @@ -76,7 +76,7 @@ function UnstructuredMesh2D(filename; RealT=Float64, periodicity=false, unsaved_ bary_weights = SVector{mesh_nnodes}(bary_weights_) arrays = (; corner_nodes, interface_info, element_node_ids, curved_check, - cornerNodeVals, tempNodes, curve_vals, + quad_corners, quad_corners_flipped, curve_values, element_is_curved, surface_curves, boundary_names) counters = (; n_corners, n_surfaces, n_elements) @@ -99,7 +99,7 @@ end function parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, cheby_nodes, bary_weights) @unpack ( corner_nodes, interface_info, element_node_ids, curved_check, - cornerNodeVals, tempNodes, curve_vals, + quad_corners, quad_corners_flipped, curve_values, element_is_curved, surface_curves, boundary_names ) = arrays @unpack n_corners, n_surfaces, n_elements = counters mesh_nnodes = length(cheby_nodes) @@ -110,7 +110,7 @@ function parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, c # readin an store the nodes that dictate the corners of the elements needed to construct the # element geometry terms for j in 1:n_corners - current_line = split(file_lines[file_idx]) + current_line = split(file_lines[file_idx]) corner_nodes[1, j] = parse(RealT, current_line[1]) corner_nodes[2, j] = parse(RealT, current_line[2]) file_idx += 1 @@ -127,7 +127,7 @@ function parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, c # container to for the interface neighbour information and connectivity n_boundaries = 0 for j in 1:n_surfaces - current_line = split(file_lines[file_idx]) + current_line = split(file_lines[file_idx]) interface_info[1, j] = parse(Int, current_line[1]) interface_info[2, j] = parse(Int, current_line[2]) interface_info[3, j] = parse(Int, current_line[3]) @@ -149,18 +149,18 @@ function parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, c for j in 1:n_elements # pull the corner node IDs - current_line = split(file_lines[file_idx]) + current_line = split(file_lines[file_idx]) element_node_ids[1, j] = parse(Int, current_line[1]) element_node_ids[2, j] = parse(Int, current_line[2]) element_node_ids[3, j] = parse(Int, current_line[3]) element_node_ids[4, j] = parse(Int, current_line[4]) for i in 1:4 # pull the (x,y) values of these corners out of the nodes array - cornerNodeVals[i, :] .= corner_nodes[:, element_node_ids[i, j]] + quad_corners[i, :] .= corner_nodes[:, element_node_ids[i, j]] end # pull the information to check if boundary is curved in order to read in additional data file_idx += 1 - current_line = split(file_lines[file_idx]) + current_line = split(file_lines[file_idx]) curved_check[1] = parse(Int, current_line[1]) curved_check[2] = parse(Int, current_line[2]) curved_check[3] = parse(Int, current_line[3]) @@ -178,32 +178,30 @@ function parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, c # flip node ordering to make sure the element is right-handed for the interpolations m1 = 1 m2 = 2 - @views tempNodes[1, :] .= cornerNodeVals[4, :] - @views tempNodes[2, :] .= cornerNodeVals[2, :] - @views tempNodes[3, :] .= cornerNodeVals[3, :] - @views tempNodes[4, :] .= cornerNodeVals[1, :] + @views quad_corners_flipped[1, :] .= quad_corners[4, :] + @views quad_corners_flipped[2, :] .= quad_corners[2, :] + @views quad_corners_flipped[3, :] .= quad_corners[3, :] + @views quad_corners_flipped[4, :] .= quad_corners[1, :] for i in 1:4 if curved_check[i] == 0 - # when curved_check[i] is 0 then the "curve" from cornerNode(i) to cornerNode(i+1) is a + # when curved_check[i] is 0 then the "curve" from corner `i` to corner `i+1` is a # straight line. So we must construct the interpolant for this line for k in 1:mesh_nnodes - curve_vals[k, 1] = tempNodes[m1, 1] + 0.5 * (cheby_nodes[k] + 1.0) * ( tempNodes[m2, 1] - - tempNodes[m1, 1]) - curve_vals[k, 2] = tempNodes[m1, 2] + 0.5 * (cheby_nodes[k] + 1.0) * ( tempNodes[m2, 2] - - tempNodes[m1, 2]) + curve_values[k, 1] = linear_interpolate(cheby_nodes[k], quad_corners_flipped[m1, 1], quad_corners_flipped[m2, 1]) + curve_values[k, 2] = linear_interpolate(cheby_nodes[k], quad_corners_flipped[m1, 2], quad_corners_flipped[m2, 2]) end else # when curved_check[i] is 1 this curved boundary information is supplied by the mesh # generator. So we just read it into a work array for k in 1:mesh_nnodes - file_idx += 1 - current_line = split(file_lines[file_idx]) - curve_vals[k, 1] = parse(RealT,current_line[1]) - curve_vals[k, 2] = parse(RealT,current_line[2]) + file_idx += 1 + current_line = split(file_lines[file_idx]) + curve_values[k, 1] = parse(RealT,current_line[1]) + curve_values[k, 2] = parse(RealT,current_line[2]) end end # construct the curve interpolant for the current side - surface_curves[i, j] = CurvedSurfaceT(cheby_nodes, bary_weights, copy(curve_vals)) + surface_curves[i, j] = CurvedSurfaceT(cheby_nodes, bary_weights, copy(curve_values)) # indexing update that contains a "flip" to ensure correct element orientation # if we need to construct the straight line "curves" when curved_check[i] == 0 m1 += 1 diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 0c9f6536f2c..b3ae29e35f1 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -261,7 +261,7 @@ function init_boundaries_iter_face_inner(info, boundaries, boundary_id, mesh) quad_id = offset + local_quad_id # Write data to boundaries container - # p4est uses zero-based indexing; convert to one-based indexing + # `p4est` uses zero-based indexing; convert to one-based indexing boundaries.neighbor_ids[boundary_id] = quad_id + 1 # Face at which the boundary lies @@ -397,7 +397,7 @@ function reinitialize_containers!(mesh::P4estMesh, equations, dg::DGSEM, cache) resize!(mortars, required.mortars) # re-initialize containers together to reduce - # the number of iterations over the mesh in p4est + # the number of iterations over the mesh in `p4est` init_surfaces!(interfaces, mortars, boundaries, mesh) end @@ -464,7 +464,7 @@ function init_surfaces_iter_face_inner(info, user_data) end function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMesh) - # Let p4est iterate over all interfaces and call init_surfaces_iter_face + # Let `p4est` iterate over all interfaces and call init_surfaces_iter_face iter_face_c = cfunction(init_surfaces_iter_face, Val(ndims(mesh))) user_data = InitSurfacesIterFaceUserData( interfaces, mortars, boundaries, mesh) @@ -492,7 +492,7 @@ function init_interfaces_iter_face_inner(info, sides, user_data) quad_ids = offsets + local_quad_ids # Write data to interfaces container - # p4est uses zero-based indexing; convert to one-based indexing + # `p4est` uses zero-based indexing; convert to one-based indexing interfaces.neighbor_ids[1, interface_id] = quad_ids[1] + 1 interfaces.neighbor_ids[2, interface_id] = quad_ids[2] + 1 @@ -526,7 +526,7 @@ function init_boundaries_iter_face_inner(info, user_data) quad_id = offset + local_quad_id # Write data to boundaries container - # p4est uses zero-based indexing; convert to one-based indexing + # `p4est` uses zero-based indexing; convert to one-based indexing boundaries.neighbor_ids[boundary_id] = quad_id + 1 # Face at which the boundary lies @@ -580,7 +580,7 @@ function init_mortars_iter_face_inner(info, sides, user_data) end # Write data to mortar container, 1 and 2 are the small elements - # p4est uses zero-based indexing; convert to one-based indexing + # `p4est` uses zero-based indexing; convert to one-based indexing mortars.neighbor_ids[1:end-1, mortar_id] .= small_quad_ids[:] .+ 1 # Last entry is the large element mortars.neighbor_ids[end, mortar_id] = large_quad_id + 1 @@ -634,7 +634,7 @@ cfunction(::typeof(count_surfaces_iter_face), ::Val{2}) = @cfunction(count_surfa cfunction(::typeof(count_surfaces_iter_face), ::Val{3}) = @cfunction(count_surfaces_iter_face, Cvoid, (Ptr{p8est_iter_face_info_t}, Ptr{Cvoid})) function count_required_surfaces(mesh::P4estMesh) - # Let p4est iterate over all interfaces and call count_surfaces_iter_face + # Let `p4est` iterate over all interfaces and call count_surfaces_iter_face iter_face_c = cfunction(count_surfaces_iter_face, Val(ndims(mesh))) # interfaces, mortars, boundaries diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index d809aa59c82..72b23f9026e 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -49,7 +49,7 @@ function calc_node_coordinates!(node_coordinates, matrix2 = similar(matrix1) baryweights_in = barycentric_weights(mesh.nodes) - # Macros from p4est + # Macros from `p4est` p4est_root_len = 1 << P4EST_MAXLEVEL p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l) diff --git a/src/solvers/dgsem_p4est/containers_3d.jl b/src/solvers/dgsem_p4est/containers_3d.jl index fc9cefe25b7..825288d73c0 100644 --- a/src/solvers/dgsem_p4est/containers_3d.jl +++ b/src/solvers/dgsem_p4est/containers_3d.jl @@ -10,7 +10,7 @@ function init_elements!(elements, mesh::P4estMesh{3}, basis::LobattoLegendreBasi @unpack node_coordinates, jacobian_matrix, contravariant_vectors, inverse_jacobian = elements - calc_node_coordinates!(node_coordinates, mesh, basis.nodes) + calc_node_coordinates!(node_coordinates, mesh, basis) for element in 1:ncells(mesh) calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) @@ -40,7 +40,7 @@ end function calc_node_coordinates!(node_coordinates, mesh::P4estMesh{3}, nodes::AbstractVector) - # Macros from p4est + # Macros from `p4est` p4est_root_len = 1 << P4EST_MAXLEVEL p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l) @@ -181,14 +181,14 @@ end end -# Convert p4est orientation code to node indices. +# Convert `p4est` orientation code to node indices. # Return node indices that index "my side" wrt "other side", # i.e., i and j are indices of other side. function orientation_to_indices_p4est(my_face, other_face, orientation_code) # my_face and other_face are the face directions (zero-based) # of "my side" and "other side" respectively. # Face corner 0 of the face with the lower face direction connects to a corner of the other face. - # The number of this corner is the orientation code in p4est. + # The number of this corner is the orientation code in `p4est`. lower = my_face <= other_face # x_pos, y_neg, and z_pos are the directions in which the face has right-handed coordinates @@ -200,7 +200,7 @@ function orientation_to_indices_p4est(my_face, other_face, orientation_code) # orientations when looked at from the same side of the interface. flipped = my_right_handed == other_right_handed - # In the folowing illustrations, p4est's face corner numbering is shown. + # In the folowing illustrations, the face corner numbering of `p4est` is shown. # ξ and η are the local coordinates of the respective face. # We're looking at both faces from the same side of the interface, so that "other side" # (in the illustrations on the left) has right-handed coordinates. diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index 516c5d45f3f..dd38fdda39d 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -9,8 +9,8 @@ # It constructs the basic `cache` used throughout the simulation to compute # the RHS etc. function create_cache(mesh::P4estMesh, equations::AbstractEquations, dg::DG, ::Any, ::Type{uEltype}) where {uEltype<:Real} - # Make sure to balance the p4est before creating any containers - # in case someone has tampered with the p4est after creating the mesh + # Make sure to balance the `p4est` before creating any containers + # in case someone has tampered with the `p4est` after creating the mesh balance!(mesh) elements = init_elements(mesh, equations, dg.basis, uEltype) diff --git a/src/solvers/dgsem_unstructured/mappings_geometry_curved_2d.jl b/src/solvers/dgsem_unstructured/mappings_geometry_curved_2d.jl index 217ce031b2f..8ad018bd08a 100644 --- a/src/solvers/dgsem_unstructured/mappings_geometry_curved_2d.jl +++ b/src/solvers/dgsem_unstructured/mappings_geometry_curved_2d.jl @@ -76,7 +76,7 @@ end # construct the (x,y) node coordinates in the volume of a curved sided element -function calc_node_coordinates!(node_coordinates, element, nodes, +function calc_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, element, nodes, surface_curves::AbstractVector{<:CurvedSurface}) for j in eachindex(nodes), i in eachindex(nodes) diff --git a/src/solvers/dgsem_unstructured/mappings_geometry_straight_2d.jl b/src/solvers/dgsem_unstructured/mappings_geometry_straight_2d.jl index ff0a0b8f7be..d4f9cf5d49d 100644 --- a/src/solvers/dgsem_unstructured/mappings_geometry_straight_2d.jl +++ b/src/solvers/dgsem_unstructured/mappings_geometry_straight_2d.jl @@ -45,7 +45,7 @@ end # construct the (x,y) node coordinates in the volume of a straight sided element -function calc_node_coordinates!(node_coordinates, element, nodes, corners) +function calc_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, element, nodes, corners) for j in eachindex(nodes), i in eachindex(nodes) node_coordinates[:, i ,j ,element] .= straight_side_quad_map(nodes[i], nodes[j], corners) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 4ce451bc4a7..5628c1bac75 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -85,6 +85,13 @@ isdir(outdir) && rm(outdir, recursive=true) tspan = (0.0, 0.3)) end + @trixi_testset "elixir_euler_wall_bc_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_wall_bc_amr.jl"), + l2 = [0.020291447969983396, 0.017479614254319948, 0.011387644425613437, 0.0514420126021293], + linf = [0.3582779022370579, 0.32073537890751663, 0.221818049107692, 0.9209559420400415], + tspan = (0.0, 0.15)) + end + @trixi_testset "elixir_eulergravity_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_convergence.jl"), l2 = [0.00024871265138964204, 0.0003370077102132591, 0.0003370077102131964, 0.0007231525513793697], diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index a3bc1364aff..94f88914e69 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -120,6 +120,12 @@ isdir(outdir) && rm(outdir, recursive=true) 0.000535219040687628], tspan = (0.0, 0.04)) end + + @trixi_testset "elixir_euler_source_terms_nonperiodic_hohqmesh.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic_hohqmesh.jl"), + l2 = [0.0042023406458005464, 0.004122532789279737, 0.0042448149597303616, 0.0036361316700401765, 0.007389845952982495], + linf = [0.04530610539892499, 0.02765695110527666, 0.05670295599308606, 0.048396544302230504, 0.1154589758186293]) + end end # Clean up afterwards: delete Trixi output directory diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index abd68a1756d..1afe0c91747 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -21,8 +21,8 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "elixir_euler_free_stream.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"), - l2 = [3.357431396258123e-14, 2.439943089578555e-13, 1.4655386790023588e-13, 4.670410488845425e-13], - linf = [1.0169198816356584e-11, 6.838458965763294e-11, 4.456759961080081e-11, 1.4207657272891083e-10], + l2 = [3.3937971107485363e-14, 2.447586447887882e-13, 1.4585205789296455e-13, 4.716993468962946e-13], + linf = [8.804734719092266e-12, 6.261270668606045e-11, 2.93670088247211e-11, 1.205400224080222e-10], tspan = (0.0, 0.1), atol = 3.0e-13) end @@ -107,8 +107,8 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "elixir_shallowwater_well_balanced.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2 = [1.2164292510839083, 2.8049631232564753e-12, 1.5664980749498454e-12, 1.216429251083908], - linf = [1.5138512282315868, 5.0263880749562876e-11, 3.5805418805003075e-11, 1.513851228231574], + l2 = [1.2164292510839076, 2.6118925543469468e-12, 1.1636046671473883e-12, 1.2164292510839079], + linf = [1.5138512282315846, 4.998482888288039e-11, 2.0246214978154587e-11, 1.513851228231574], tspan = (0.0, 0.25)) end @@ -121,8 +121,8 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "elixir_shallowwater_dirichlet.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_dirichlet.jl"), - l2 = [1.1577518608963063e-5, 4.371173779970303e-13, 4.2152984234036224e-13, 1.1577518608935235e-5], - linf = [8.394063878491842e-5, 9.632644747754261e-11, 9.54898901893391e-11, 8.394063879602065e-5], + l2 = [1.1577518608940115e-5, 4.867189932537344e-13, 4.647273240470541e-13, 1.1577518608933468e-5], + linf = [8.394063878602864e-5, 1.1469760027632646e-10, 1.1146619484429974e-10, 8.394063879602065e-5], tspan = (0.0, 2.0)) end