From 02f40d33347d657cd964585aa6ed9fb08c856e3c Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Fri, 2 Feb 2024 11:58:34 +0100 Subject: [PATCH] format (#54) --- .JuliaFormatter.toml | 8 + .github/workflows/FormatCheck.yml | 38 ++ examples/t8_step1_coarsemesh.jl | 28 +- examples/t8_step2_uniform_forest.jl | 22 +- examples/t8_step3_common.jl | 117 +++-- examples/t8_step4_partition_balance_ghost.jl | 110 +++-- examples/t8_step5_element_data.jl | 272 ++++++----- examples/t8_step6_stencil.jl | 489 +++++++++---------- examples/t8_tutorial_build_cmesh.jl | 470 +++++++++--------- src/T8code.jl | 18 +- test/cmesh/test_readmshfile.jl | 210 ++++---- test/forest/test_element_volume.jl | 95 ++-- test/runtests.jl | 32 +- test/test_all.jl | 14 +- test/test_init.jl | 12 +- test/test_refcount.jl | 106 ++-- 16 files changed, 1035 insertions(+), 1006 deletions(-) create mode 100644 .JuliaFormatter.toml create mode 100644 .github/workflows/FormatCheck.yml diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 0000000..8518d20 --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1,8 @@ +# Use SciML style: https://github.com/SciML/SciMLStyle +style = "sciml" + +# Python style alignment. See https://github.com/domluna/JuliaFormatter.jl/pull/732. +yas_style_nesting = true + +# Align struct fields for better readability of large struct definitions +align_struct_field = true diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml new file mode 100644 index 0000000..4d0dd92 --- /dev/null +++ b/.github/workflows/FormatCheck.yml @@ -0,0 +1,38 @@ +name: format-check + +on: + push: + branches: + - 'main' + tags: '*' + pull_request: + +jobs: + check-format: + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: [1] + julia-arch: [x86] + os: [ubuntu-latest] + steps: + - uses: julia-actions/setup-julia@latest + with: + version: ${{ matrix.julia-version }} + - uses: actions/checkout@v4 + - uses: julia-actions/cache@v1 + - name: Install JuliaFormatter and format + run: | + julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter"))' + julia -e 'using JuliaFormatter; format(["examples", "src/T8code.jl", "test"])' + - name: Format check + run: | + julia -e ' + out = Cmd(`git diff --name-only`) |> read |> String + if out == "" + exit(0) + else + @error "Some files have not been formatted !!!" + write(stdout, out) + exit(1) + end' diff --git a/examples/t8_step1_coarsemesh.jl b/examples/t8_step1_coarsemesh.jl index cff94db..0e252f9 100644 --- a/examples/t8_step1_coarsemesh.jl +++ b/examples/t8_step1_coarsemesh.jl @@ -23,20 +23,20 @@ using T8code.Libt8: SC_LP_PRODUCTION # \param [in] comm MPI Communicator to use. # \return The coarse mesh. function t8_step1_build_tetcube_coarse_mesh(comm) - # Build a coarse mesh of 6 tetrahedral trees that form a cube. - # You can modify the first parameter to build a cube with different - # tree shapes, i.e. T8_ECLASS_QUAD for a unit square with 1 quadrilateral tree. - # See t8_eclass.h, t8_cmesh.h for all possible shapes. - # - # The second argument is the MPI communicator to use for this cmesh. - # The remaining arguments are 3 flags that control - # do_bcast - If non-zero only the root process will build the cmesh and will broadcast it to the other processes. The result is the same. - # do_partition - If non-zero the cmesh will be partitioned among the processes. If 0 each process has a copy of the whole cmesh. - # periodic - If non-zero the cube will have periodic boundaries. That is, i.e. the left face is connected to the right face. + # Build a coarse mesh of 6 tetrahedral trees that form a cube. + # You can modify the first parameter to build a cube with different + # tree shapes, i.e. T8_ECLASS_QUAD for a unit square with 1 quadrilateral tree. + # See t8_eclass.h, t8_cmesh.h for all possible shapes. + # + # The second argument is the MPI communicator to use for this cmesh. + # The remaining arguments are 3 flags that control + # do_bcast - If non-zero only the root process will build the cmesh and will broadcast it to the other processes. The result is the same. + # do_partition - If non-zero the cmesh will be partitioned among the processes. If 0 each process has a copy of the whole cmesh. + # periodic - If non-zero the cube will have periodic boundaries. That is, i.e. the left face is connected to the right face. - cmesh = t8_cmesh_new_hypercube(T8_ECLASS_TET, comm, 0, 0, 0) + cmesh = t8_cmesh_new_hypercube(T8_ECLASS_TET, comm, 0, 0, 0) - return cmesh + return cmesh end # Write vtk (or more accurately vtu) files of the cmesh. @@ -47,13 +47,13 @@ end # If the coarse mesh would be repartitioned, then it would write the .pvtu file # and additionally one file prefix_MPIRANK.vtu per MPI rank. function t8_step1_write_cmesh_vtk(cmesh, prefix) - t8_cmesh_vtk_write_file(cmesh, prefix, 1.0) + t8_cmesh_vtk_write_file(cmesh, prefix, 1.0) end # Destroy a cmesh. This will free all allocated memory. # \param [in] cmesh A cmesh. function t8_step1_destroy_cmesh(cmesh) - t8_cmesh_destroy(Ref(cmesh)) + t8_cmesh_destroy(Ref(cmesh)) end # The prefix for our output files. diff --git a/examples/t8_step2_uniform_forest.jl b/examples/t8_step2_uniform_forest.jl index 75e278c..6be5fd6 100644 --- a/examples/t8_step2_uniform_forest.jl +++ b/examples/t8_step2_uniform_forest.jl @@ -54,12 +54,12 @@ using T8code.Libt8: SC_LP_PRODUCTION # \param [in] comm MPI Communicator to use. # \return The coarse mesh. function t8_step2_build_prismcube_coarse_mesh(comm) - # Build a coarse mesh of 2 prism trees that form a cube. - cmesh = t8_cmesh_new_hypercube(T8_ECLASS_PRISM, comm, 0, 0, 0) + # Build a coarse mesh of 2 prism trees that form a cube. + cmesh = t8_cmesh_new_hypercube(T8_ECLASS_PRISM, comm, 0, 0, 0) - t8_global_productionf(" [step2] Constructed coarse mesh with 2 prism trees.\n") + t8_global_productionf(" [step2] Constructed coarse mesh with 2 prism trees.\n") - return cmesh + return cmesh end # Build a uniform forest on a cmesh @@ -70,12 +70,12 @@ end # \return A uniform forest with the given refinement level that is # partitioned across the processes in \a comm. function t8_step2_build_uniform_forest(comm, cmesh, level) - # /* Create the refinement scheme. */ - scheme = t8_scheme_new_default_cxx() - # /* Creat the uniform forest. */ - forest = t8_forest_new_uniform(cmesh, scheme, level, 0, comm) + # /* Create the refinement scheme. */ + scheme = t8_scheme_new_default_cxx() + # /* Creat the uniform forest. */ + forest = t8_forest_new_uniform(cmesh, scheme, level, 0, comm) - return forest + return forest end # Write vtk (or more accurately vtu) files of the forest. @@ -85,7 +85,7 @@ end # This will create the file prefix.pvtu # and additionally one file prefix_MPIRANK.vtu per MPI rank. function t8_step2_write_forest_vtk(forest, prefix) - t8_forest_write_vtk(forest, prefix) + t8_forest_write_vtk(forest, prefix) end # Destroy a forest. This will free all allocated memory. @@ -95,7 +95,7 @@ end # If we do not want this behaviour, but want to reuse for example the cmesh, # we need to call t8_cmesh_ref (cmesh) before passing it to t8_forest_new_uniform. function t8_step2_destroy_forest(forest) - t8_forest_unref(Ref(forest)) + t8_forest_unref(Ref(forest)) end # The prefix for our output files. diff --git a/examples/t8_step3_common.jl b/examples/t8_step3_common.jl index c79743a..bc263b6 100644 --- a/examples/t8_step3_common.jl +++ b/examples/t8_step3_common.jl @@ -1,9 +1,9 @@ # This is our own defined data that we will pass on to the # adaptation callback. mutable struct t8_step3_adapt_data_t - midpoint :: NTuple{3,Cdouble} - refine_if_inside_radius :: Cdouble - coarsen_if_outside_radius :: Cdouble + midpoint :: NTuple{3, Cdouble} + refine_if_inside_radius :: Cdouble + coarsen_if_outside_radius :: Cdouble end # The adaptation callback function. This function will be called once for each element @@ -26,81 +26,80 @@ end # \param [in] num_elements The number of entries in \a elements elements that are defined. # \param [in] elements The element or family of elements to consider for refinement/coarsening. function t8_step3_adapt_callback(forest, forest_from, which_tree, lelement_id, - ts, is_family, num_elements, elements_ptr) :: Cint - # Our adaptation criterion is to look at the midpoint coordinates of the current element and if - # they are inside a sphere around a given midpoint we refine, if they are outside, we coarsen. + ts, is_family, num_elements, elements_ptr)::Cint + # Our adaptation criterion is to look at the midpoint coordinates of the current element and if + # they are inside a sphere around a given midpoint we refine, if they are outside, we coarsen. - centroid = Vector{Cdouble}(undef,3) # Will hold the element midpoint. - # In t8_step3_adapt_forest we pass a t8_step3_adapt_data pointer as user data to the - # t8_forest_new_adapt function. This pointer is stored as the used data of the new forest - # and we can now access it with t8_forest_get_user_data (forest). - adapt_data_ptr = Ptr{t8_step3_adapt_data_t}(t8_forest_get_user_data(forest)) + centroid = Vector{Cdouble}(undef, 3) # Will hold the element midpoint. + # In t8_step3_adapt_forest we pass a t8_step3_adapt_data pointer as user data to the + # t8_forest_new_adapt function. This pointer is stored as the used data of the new forest + # and we can now access it with t8_forest_get_user_data (forest). + adapt_data_ptr = Ptr{t8_step3_adapt_data_t}(t8_forest_get_user_data(forest)) - # You can use assert for assertions that are active in debug mode (when configured with --enable-debug). - # If the condition is not true, then the code will abort. - # In this case, we want to make sure that we actually did set a user pointer to forest and thus - # did not get the NULL pointer from t8_forest_get_user_data. - @T8_ASSERT(adapt_data_ptr != C_NULL) + # You can use assert for assertions that are active in debug mode (when configured with --enable-debug). + # If the condition is not true, then the code will abort. + # In this case, we want to make sure that we actually did set a user pointer to forest and thus + # did not get the NULL pointer from t8_forest_get_user_data. + @T8_ASSERT(adapt_data_ptr!=C_NULL) - adapt_data = unsafe_load(adapt_data_ptr) + adapt_data = unsafe_load(adapt_data_ptr) - elements = unsafe_wrap(Array, elements_ptr, num_elements) + elements = unsafe_wrap(Array, elements_ptr, num_elements) - # Compute the element's centroid coordinates. - t8_forest_element_centroid(forest_from, which_tree, elements[1], pointer(centroid)) + # Compute the element's centroid coordinates. + t8_forest_element_centroid(forest_from, which_tree, elements[1], pointer(centroid)) - # Compute the distance to our sphere midpoint. - dist = sqrt(sum((centroid .- adapt_data.midpoint).^2)) - if dist < adapt_data.refine_if_inside_radius - # Refine this element. - return 1 - elseif is_family == 1 && dist > adapt_data.coarsen_if_outside_radius - # Coarsen this family. Note that we check for is_family before, since returning < 0 - # if we do not have a family as input is illegal. - return -1 - end + # Compute the distance to our sphere midpoint. + dist = sqrt(sum((centroid .- adapt_data.midpoint) .^ 2)) + if dist < adapt_data.refine_if_inside_radius + # Refine this element. + return 1 + elseif is_family == 1 && dist > adapt_data.coarsen_if_outside_radius + # Coarsen this family. Note that we check for is_family before, since returning < 0 + # if we do not have a family as input is illegal. + return -1 + end - # Do not change this element. - return 0 + # Do not change this element. + return 0 end # Adapt a forest according to our t8_step3_adapt_callback function. # This will create a new forest and return it. function t8_step3_adapt_forest(forest) - adapt_data = t8_step3_adapt_data_t( - (0.5, 0.5, 1.0), # Midpoints of the sphere. - 0.2, # Refine if inside this radius. - 0.4 # Coarsen if outside this radius. - ) + adapt_data = t8_step3_adapt_data_t((0.5, 0.5, 1.0), # Midpoints of the sphere. + 0.2, # Refine if inside this radius. + 0.4) - # Check that forest is a committed, that is valid and usable, forest. - @T8_ASSERT(t8_forest_is_committed(forest) == 1) + # Check that forest is a committed, that is valid and usable, forest. + @T8_ASSERT(t8_forest_is_committed(forest)==1) - # Create a new forest that is adapted from \a forest with our adaptation callback. - # We provide the adapt_data as user data that is stored as the used_data pointer of the - # new forest (see also t8_forest_set_user_data). - # The 0, 0 arguments are flags that control - # recursive - If non-zero adaptation is recursive, thus if an element is adapted the children - # or parents are plugged into the callback again recursively until the forest does not - # change any more. If you use this you should ensure that refinement will stop eventually. - # One way is to check the element's level against a given maximum level. - # do_face_ghost - If non-zero additionally a layer of ghost elements is created for the forest. - # We will discuss ghost in later steps of the tutorial. - forest_adapt = t8_forest_new_adapt(forest, @t8_adapt_callback(t8_step3_adapt_callback), 0, 0, Ref(adapt_data)) + # Create a new forest that is adapted from \a forest with our adaptation callback. + # We provide the adapt_data as user data that is stored as the used_data pointer of the + # new forest (see also t8_forest_set_user_data). + # The 0, 0 arguments are flags that control + # recursive - If non-zero adaptation is recursive, thus if an element is adapted the children + # or parents are plugged into the callback again recursively until the forest does not + # change any more. If you use this you should ensure that refinement will stop eventually. + # One way is to check the element's level against a given maximum level. + # do_face_ghost - If non-zero additionally a layer of ghost elements is created for the forest. + # We will discuss ghost in later steps of the tutorial. + forest_adapt = t8_forest_new_adapt(forest, @t8_adapt_callback(t8_step3_adapt_callback), + 0, 0, Ref(adapt_data)) - return forest_adapt + return forest_adapt end # Print the local and global number of elements of a forest. function t8_step3_print_forest_information(forest) - # Check that forest is a committed, that is valid and usable, forest. - @T8_ASSERT(t8_forest_is_committed(forest) == 1) + # Check that forest is a committed, that is valid and usable, forest. + @T8_ASSERT(t8_forest_is_committed(forest)==1) - # Get the local number of elements. - local_num_elements = t8_forest_get_local_num_elements(forest) - # Get the global number of elements. - global_num_elements = t8_forest_get_global_num_elements(forest) + # Get the local number of elements. + local_num_elements = t8_forest_get_local_num_elements(forest) + # Get the global number of elements. + global_num_elements = t8_forest_get_global_num_elements(forest) - t8_global_productionf(" [step3] Local number of elements:\t\t%i\n", local_num_elements) - t8_global_productionf(" [step3] Global number of elements:\t%li\n", global_num_elements) + t8_global_productionf(" [step3] Local number of elements:\t\t%i\n", local_num_elements) + t8_global_productionf(" [step3] Global number of elements:\t%li\n", global_num_elements) end diff --git a/examples/t8_step4_partition_balance_ghost.jl b/examples/t8_step4_partition_balance_ghost.jl index a942f31..9966f1d 100644 --- a/examples/t8_step4_partition_balance_ghost.jl +++ b/examples/t8_step4_partition_balance_ghost.jl @@ -113,65 +113,66 @@ using T8code.Libt8: SC_LP_PRODUCTION # In this function we create a new forest that repartitions a given forest # and has a layer of ghost elements. function t8_step4_partition_ghost(forest) - # Check that forest is a committed, that is a valid and usable, forest. - @T8_ASSERT(t8_forest_is_committed(forest) == 1) - - # Initialize. - new_forest_ref = Ref(t8_forest_t()) - t8_forest_init(new_forest_ref) - new_forest = new_forest_ref[] - - # Tell the new_forest that is should partition the existing forest. - # This will change the distribution of the forest elements among the processes - # in such a way that afterwards each process has the same number of elements - # (+- 1 if the number of elements is not divisible by the number of processes). - # - # The third 0 argument is the flag 'partition_for_coarsening' which is currently not - # implemented. Once it is, this will ensure that a family of elements will not be split - # across multiple processes and thus one level coarsening is always possible (see also the - # comments on coarsening in t8_step3). - t8_forest_set_partition(new_forest, forest, 0) - - # Tell the new_forest to create a ghost layer. - # This will gather those face neighbor elements of process local element that reside - # on a different process. - # - # We currently support ghost mode T8_GHOST_FACES that creates face neighbor ghost elements - # and will in future also support other modes for edge/vertex neighbor ghost elements. - t8_forest_set_ghost(new_forest, 1, T8_GHOST_FACES) - - # Commit the forest, this step will perform the partitioning and ghost layer creation. - t8_forest_commit(new_forest) - - return new_forest + # Check that forest is a committed, that is a valid and usable, forest. + @T8_ASSERT(t8_forest_is_committed(forest)==1) + + # Initialize. + new_forest_ref = Ref(t8_forest_t()) + t8_forest_init(new_forest_ref) + new_forest = new_forest_ref[] + + # Tell the new_forest that is should partition the existing forest. + # This will change the distribution of the forest elements among the processes + # in such a way that afterwards each process has the same number of elements + # (+- 1 if the number of elements is not divisible by the number of processes). + # + # The third 0 argument is the flag 'partition_for_coarsening' which is currently not + # implemented. Once it is, this will ensure that a family of elements will not be split + # across multiple processes and thus one level coarsening is always possible (see also the + # comments on coarsening in t8_step3). + t8_forest_set_partition(new_forest, forest, 0) + + # Tell the new_forest to create a ghost layer. + # This will gather those face neighbor elements of process local element that reside + # on a different process. + # + # We currently support ghost mode T8_GHOST_FACES that creates face neighbor ghost elements + # and will in future also support other modes for edge/vertex neighbor ghost elements. + t8_forest_set_ghost(new_forest, 1, T8_GHOST_FACES) + + # Commit the forest, this step will perform the partitioning and ghost layer creation. + t8_forest_commit(new_forest) + + return new_forest end # In this function we adapt a forest as in step3 and balance it. In our main # program the input forest is already adapted and then the resulting twice # adapted forest will be unbalanced. function t8_step4_balance(forest) - # Adapt the input forest. - unbalanced_forest = t8_step3_adapt_forest(forest) - - # Output to vtk. - t8_forest_write_vtk(unbalanced_forest, "t8_step4_unbalanced_forest") - t8_global_productionf(" [step4] Wrote unbalanced forest to vtu files: %s*\n", "t8_step4_unbalanced_forest") - - # Initialize new forest. - balanced_forest_ref = Ref(t8_forest_t()) - t8_forest_init(balanced_forest_ref) - balanced_forest = balanced_forest_ref[] - - # Specify that this forest should result from balancing unbalanced_forest. - # The last argument is the flag 'no_repartition'. - # Since balancing will refine elements, the load-balance will be broken afterwards. - # Setting this flag to false (no_repartition = false -> yes repartition) will repartition - # the forest after balance, such that every process has the same number of elements afterwards. - t8_forest_set_balance(balanced_forest, unbalanced_forest, 0) - # Commit the forest. - t8_forest_commit(balanced_forest) - - return balanced_forest + # Adapt the input forest. + unbalanced_forest = t8_step3_adapt_forest(forest) + + # Output to vtk. + t8_forest_write_vtk(unbalanced_forest, "t8_step4_unbalanced_forest") + t8_global_productionf(" [step4] Wrote unbalanced forest to vtu files: %s*\n", + "t8_step4_unbalanced_forest") + + # Initialize new forest. + balanced_forest_ref = Ref(t8_forest_t()) + t8_forest_init(balanced_forest_ref) + balanced_forest = balanced_forest_ref[] + + # Specify that this forest should result from balancing unbalanced_forest. + # The last argument is the flag 'no_repartition'. + # Since balancing will refine elements, the load-balance will be broken afterwards. + # Setting this flag to false (no_repartition = false -> yes repartition) will repartition + # the forest after balance, such that every process has the same number of elements afterwards. + t8_forest_set_balance(balanced_forest, unbalanced_forest, 0) + # Commit the forest. + t8_forest_commit(balanced_forest) + + return balanced_forest end include("t8_step3_common.jl") @@ -221,7 +222,8 @@ t8_step3_print_forest_information(forest); # Write forest to vtu files. t8_forest_write_vtk(forest, prefix_uniform) -t8_global_productionf(" [step4] Wrote uniform level %i forest to vtu files: %s*\n", level, prefix_uniform) +t8_global_productionf(" [step4] Wrote uniform level %i forest to vtu files: %s*\n", level, + prefix_uniform) # # Adapt the forest. diff --git a/examples/t8_step5_element_data.jl b/examples/t8_step5_element_data.jl index b094967..1130287 100644 --- a/examples/t8_step5_element_data.jl +++ b/examples/t8_step5_element_data.jl @@ -53,102 +53,102 @@ include("t8_step3_common.jl") # The data that we want to store for each element. # In this example we want to store the element's level and volume. */ struct t8_step5_data_per_element_t - level :: Cint - volume :: Cdouble + level :: Cint + volume :: Cdouble end function t8_step5_build_forest(comm, level) - cmesh = t8_cmesh_new_hypercube_hybrid(comm, 0, 0) - scheme = t8_scheme_new_default_cxx() - - adapt_data = t8_step3_adapt_data_t( - (0.5, 0.5, 1.0), # Midpoints of the sphere. - 0.2, # Refine if inside this radius. - 0.4 # Coarsen if outside this radius. - ) - - # Start with a uniform forest. - forest = t8_forest_new_uniform(cmesh, scheme, level, 0, comm) - - # Adapt, partition, balance and create ghost elements all in the same step. - forest_apbg_ref = Ref(t8_forest_t()) - t8_forest_init(forest_apbg_ref) - forest_apbg = forest_apbg_ref[] - - t8_forest_set_user_data(forest_apbg, Ref(adapt_data)) - t8_forest_set_adapt(forest_apbg, forest, @t8_adapt_callback(t8_step3_adapt_callback), 0) - t8_forest_set_partition(forest_apbg, C_NULL, 0) - t8_forest_set_balance(forest_apbg, C_NULL, 0) - t8_forest_set_ghost(forest_apbg, 1, T8_GHOST_FACES) - t8_forest_commit(forest_apbg) - - return forest_apbg + cmesh = t8_cmesh_new_hypercube_hybrid(comm, 0, 0) + scheme = t8_scheme_new_default_cxx() + + adapt_data = t8_step3_adapt_data_t((0.5, 0.5, 1.0), # Midpoints of the sphere. + 0.2, # Refine if inside this radius. + 0.4) + + # Start with a uniform forest. + forest = t8_forest_new_uniform(cmesh, scheme, level, 0, comm) + + # Adapt, partition, balance and create ghost elements all in the same step. + forest_apbg_ref = Ref(t8_forest_t()) + t8_forest_init(forest_apbg_ref) + forest_apbg = forest_apbg_ref[] + + t8_forest_set_user_data(forest_apbg, Ref(adapt_data)) + t8_forest_set_adapt(forest_apbg, forest, @t8_adapt_callback(t8_step3_adapt_callback), 0) + t8_forest_set_partition(forest_apbg, C_NULL, 0) + t8_forest_set_balance(forest_apbg, C_NULL, 0) + t8_forest_set_ghost(forest_apbg, 1, T8_GHOST_FACES) + t8_forest_commit(forest_apbg) + + return forest_apbg end function t8_step5_create_element_data(forest) - # Check that forest is a committed, that is valid and usable, forest. - @T8_ASSERT(t8_forest_is_committed(forest) == 1) - - # Get the number of local elements of forest. - num_local_elements = t8_forest_get_local_num_elements(forest) - # Get the number of ghost elements of forest. - num_ghost_elements = t8_forest_get_num_ghosts(forest) - - # Now we need to build an array of our data that is as long as the number of - # elements plus the number of ghosts. You can use any allocator such as new, - # malloc or the t8code provide allocation macro T8_ALLOC. Note that in the - # latter case you need to use T8_FREE in order to free the memory. - element_data = Array{t8_step5_data_per_element_t}(undef, num_local_elements + num_ghost_elements) - - # Note: We will later need to associate this data with an sc_array in order to exchange the values for - # the ghost elements, which we can do with sc_array_new_data (see t8_step5_exchange_ghost_data). - # We could also have directly allocated the data here in an sc_array with - # sc_array_new_count (sizeof (struct data_per_element), num_local_elements + num_ghost_elements) - - # Let us now fill the data with something. For this, we iterate through all - # trees and for each tree through all its elements, calling - # t8_forest_get_element_in_tree to get a pointer to the current element. - # This is the recommended and most performant way. An alternative is to - # iterate over the number of local elements and use t8_forest_get_element. - # However, this function needs to perform a binary search for the element and - # the tree it is in, while t8_forest_get_element_in_tree has a constant look - # up time. You should only use t8_forest_get_element if you do not know in - # which tree an element is. - - # Get the number of trees that have elements of this process. - num_local_trees = t8_forest_get_num_local_trees(forest) - - current_index = 0 - for itree = 0:num_local_trees-1 - # This loop iterates through all local trees in the forest. - # Each tree may have a different element class (quad/tri/hex/tet etc.) and therefore - # also a different way to interpret its elements. In order to be able to handle elements - # of a tree, we need to get its eclass_scheme, and in order to so we first get its eclass. - tree_class = t8_forest_get_tree_class(forest, itree) - eclass_scheme = t8_forest_get_eclass_scheme(forest, tree_class) - - # Get the number of elements of this tree. - num_elements_in_tree = t8_forest_get_tree_num_elements(forest, itree) - # This loop iterates through all the local elements of the forest in the current tree. - for ielement = 0:num_elements_in_tree-1 - current_index += 1 # Note: Julia has 1-based indexing, while C/C++ starts with 0. - - # We can now write to the position current_index into our array in order - # to store data for this element. */ Since in this example we want to - # compute the data based on the element in question, we need to get a - # pointer to this element. - element = t8_forest_get_element_in_tree(forest, itree, ielement) - - # We want to store the elements level and its volume as data. We compute these - # via the eclass_scheme and the forest_element interface. - level = t8_element_level(eclass_scheme, element) - volume = t8_forest_element_volume(forest, itree, element) - - element_data[current_index] = t8_step5_data_per_element_t(level,volume) + # Check that forest is a committed, that is valid and usable, forest. + @T8_ASSERT(t8_forest_is_committed(forest)==1) + + # Get the number of local elements of forest. + num_local_elements = t8_forest_get_local_num_elements(forest) + # Get the number of ghost elements of forest. + num_ghost_elements = t8_forest_get_num_ghosts(forest) + + # Now we need to build an array of our data that is as long as the number of + # elements plus the number of ghosts. You can use any allocator such as new, + # malloc or the t8code provide allocation macro T8_ALLOC. Note that in the + # latter case you need to use T8_FREE in order to free the memory. + element_data = Array{t8_step5_data_per_element_t}(undef, + num_local_elements + + num_ghost_elements) + + # Note: We will later need to associate this data with an sc_array in order to exchange the values for + # the ghost elements, which we can do with sc_array_new_data (see t8_step5_exchange_ghost_data). + # We could also have directly allocated the data here in an sc_array with + # sc_array_new_count (sizeof (struct data_per_element), num_local_elements + num_ghost_elements) + + # Let us now fill the data with something. For this, we iterate through all + # trees and for each tree through all its elements, calling + # t8_forest_get_element_in_tree to get a pointer to the current element. + # This is the recommended and most performant way. An alternative is to + # iterate over the number of local elements and use t8_forest_get_element. + # However, this function needs to perform a binary search for the element and + # the tree it is in, while t8_forest_get_element_in_tree has a constant look + # up time. You should only use t8_forest_get_element if you do not know in + # which tree an element is. + + # Get the number of trees that have elements of this process. + num_local_trees = t8_forest_get_num_local_trees(forest) + + current_index = 0 + for itree in 0:(num_local_trees - 1) + # This loop iterates through all local trees in the forest. + # Each tree may have a different element class (quad/tri/hex/tet etc.) and therefore + # also a different way to interpret its elements. In order to be able to handle elements + # of a tree, we need to get its eclass_scheme, and in order to so we first get its eclass. + tree_class = t8_forest_get_tree_class(forest, itree) + eclass_scheme = t8_forest_get_eclass_scheme(forest, tree_class) + + # Get the number of elements of this tree. + num_elements_in_tree = t8_forest_get_tree_num_elements(forest, itree) + # This loop iterates through all the local elements of the forest in the current tree. + for ielement in 0:(num_elements_in_tree - 1) + current_index += 1 # Note: Julia has 1-based indexing, while C/C++ starts with 0. + + # We can now write to the position current_index into our array in order + # to store data for this element. */ Since in this example we want to + # compute the data based on the element in question, we need to get a + # pointer to this element. + element = t8_forest_get_element_in_tree(forest, itree, ielement) + + # We want to store the elements level and its volume as data. We compute these + # via the eclass_scheme and the forest_element interface. + level = t8_element_level(eclass_scheme, element) + volume = t8_forest_element_volume(forest, itree, element) + + element_data[current_index] = t8_step5_data_per_element_t(level, volume) + end end - end - return element_data + return element_data end # Each process has computed the data entries for its local elements. In order @@ -157,15 +157,17 @@ end # entries of our element data array with the value on the process that owns the # corresponding element. */ function t8_step5_exchange_ghost_data(forest, element_data) - # t8_forest_ghost_exchange_data expects an sc_array (of length num_local_elements + num_ghosts). - # We wrap our data array to an sc_array. - sc_array_wrapper = sc_array_new_data(pointer(element_data), sizeof(t8_step5_data_per_element_t), length(element_data)) + # t8_forest_ghost_exchange_data expects an sc_array (of length num_local_elements + num_ghosts). + # We wrap our data array to an sc_array. + sc_array_wrapper = sc_array_new_data(pointer(element_data), + sizeof(t8_step5_data_per_element_t), + length(element_data)) - # Carry out the data exchange. The entries with indices > num_local_elements will get overwritten. - t8_forest_ghost_exchange_data(forest, sc_array_wrapper) + # Carry out the data exchange. The entries with indices > num_local_elements will get overwritten. + t8_forest_ghost_exchange_data(forest, sc_array_wrapper) - # Destroy the wrapper array. This will not free the data memory since we used sc_array_new_data. - sc_array_destroy(sc_array_wrapper) + # Destroy the wrapper array. This will not free the data memory since we used sc_array_new_data. + sc_array_destroy(sc_array_wrapper) end # Write the forest as vtu and also write the element's volumes in the file. @@ -176,38 +178,36 @@ end # We support two types: T8_VTK_SCALAR - One double per element # and T8_VTK_VECTOR - 3 doubles per element function t8_step5_output_data_to_vtu(forest, element_data, prefix) - num_elements = t8_forest_get_local_num_elements(forest) - # We need to allocate a new array to store the volumes on their own. - # This array has one entry per local element. */ - element_volumes = Vector{Cdouble}(undef, num_elements) - - # Copy the element's volumes from our data array to the output array. - for ielem = 1:num_elements - element_volumes[ielem] = element_data[ielem].volume - end - - # The number of user defined data fields to write. - num_data = 1 - - # WARNING: This code hangs for Julia v1.8.* or older. Use at least Julia v1.9. - # For each user defined data field we need one t8_vtk_data_field_t variable. - vtk_data = t8_vtk_data_field_t( - T8_VTK_SCALAR, # Set the type of this variable. Since we have one value per element, we pick T8_VTK_SCALAR. - NTuple{8192, Cchar}(rpad("Element volume\0", 8192, ' ')), # The name of the field as should be written to the file. - pointer(element_volumes), # Pointer to the data. - ) - - # To write user defined data, we need to extended output function - # t8_forest_vtk_write_file from t8_forest_vtk.h. Despite writin user data, - # it also offers more control over which properties of the forest to write. - write_treeid = 1 - write_mpirank = 1 - write_level = 1 - write_element_id = 1 - write_ghosts = 0 - t8_forest_write_vtk_ext(forest, prefix, write_treeid, write_mpirank, - write_level, write_element_id, write_ghosts, - 0, 0, num_data, Ref(vtk_data)) + num_elements = t8_forest_get_local_num_elements(forest) + # We need to allocate a new array to store the volumes on their own. + # This array has one entry per local element. */ + element_volumes = Vector{Cdouble}(undef, num_elements) + + # Copy the element's volumes from our data array to the output array. + for ielem in 1:num_elements + element_volumes[ielem] = element_data[ielem].volume + end + + # The number of user defined data fields to write. + num_data = 1 + + # WARNING: This code hangs for Julia v1.8.* or older. Use at least Julia v1.9. + # For each user defined data field we need one t8_vtk_data_field_t variable. + vtk_data = t8_vtk_data_field_t(T8_VTK_SCALAR, # Set the type of this variable. Since we have one value per element, we pick T8_VTK_SCALAR. + NTuple{8192, Cchar}(rpad("Element volume\0", 8192, ' ')), # The name of the field as should be written to the file. + pointer(element_volumes)) + + # To write user defined data, we need to extended output function + # t8_forest_vtk_write_file from t8_forest_vtk.h. Despite writin user data, + # it also offers more control over which properties of the forest to write. + write_treeid = 1 + write_mpirank = 1 + write_level = 1 + write_element_id = 1 + write_ghosts = 0 + t8_forest_write_vtk_ext(forest, prefix, write_treeid, write_mpirank, + write_level, write_element_id, write_ghosts, + 0, 0, num_data, Ref(vtk_data)) end # The prefix for our output files. @@ -256,8 +256,9 @@ element_data = t8_step5_create_element_data(forest) t8_global_productionf(" [step5] Computed level and volume data for local elements.\n") if t8_forest_get_local_num_elements(forest) > 0 - # Output the stored data of the first local element (if it exists). - t8_global_productionf(" [step5] Element 0 has level %i and volume %e.\n", element_data[1].level, element_data[1].volume) + # Output the stored data of the first local element (if it exists). + t8_global_productionf(" [step5] Element 0 has level %i and volume %e.\n", + element_data[1].level, element_data[1].volume) end # @@ -267,16 +268,17 @@ t8_step5_exchange_ghost_data(forest, element_data) t8_global_productionf(" [step5] Exchanged ghost data.\n") if t8_forest_get_num_ghosts(forest) > 0 - # Output the data of the first ghost element (if it exists). - first_ghost_index = t8_forest_get_local_num_elements(forest) - t8_global_productionf(" [step5] Ghost 0 has level %i and volume %e.\n", - element_data[first_ghost_index + 1].level, - element_data[first_ghost_index + 1].volume) + # Output the data of the first ghost element (if it exists). + first_ghost_index = t8_forest_get_local_num_elements(forest) + t8_global_productionf(" [step5] Ghost 0 has level %i and volume %e.\n", + element_data[first_ghost_index + 1].level, + element_data[first_ghost_index + 1].volume) end # Output the volume data to vtu. t8_step5_output_data_to_vtu(forest, element_data, prefix_forest_with_data) -t8_global_productionf(" [step5] Wrote forest and volume data to %s*.\n", prefix_forest_with_data) +t8_global_productionf(" [step5] Wrote forest and volume data to %s*.\n", + prefix_forest_with_data) # # Clean-up. diff --git a/examples/t8_step6_stencil.jl b/examples/t8_step6_stencil.jl index dd12c9e..265a62d 100644 --- a/examples/t8_step6_stencil.jl +++ b/examples/t8_step6_stencil.jl @@ -49,25 +49,25 @@ include("t8_step3_common.jl") # The data that we want to store for each element. struct t8_step6_data_per_element_t - # The first three data fields are not necessary for our - # computations in this step, but are left here for reference. + # The first three data fields are not necessary for our + # computations in this step, but are left here for reference. - level :: Cint + level::Cint - volume :: Cdouble - midpoint :: NTuple{3,Cdouble} + volume::Cdouble + midpoint::NTuple{3, Cdouble} - # Element length in x- and y-direction. - dx :: Cdouble - dy :: Cdouble + # Element length in x- and y-direction. + dx::Cdouble + dy::Cdouble - # `Height` which is filled according to the position of the element. - # in the computational domain. - height :: Cdouble + # `Height` which is filled according to the position of the element. + # in the computational domain. + height::Cdouble - # Storage for our finite difference computations. - schlieren :: Cdouble - curvature :: Cdouble + # Storage for our finite difference computations. + schlieren::Cdouble + curvature::Cdouble end # In this function we first allocate a new uniformly refined forest at given @@ -76,221 +76,224 @@ end # properties of the first ("root") forest and deallocates it. The final # adapted and committed forest is returned back to the calling scope. function t8_step6_build_forest(comm, dim, level) - cmesh = t8_cmesh_new_periodic(comm, dim) + cmesh = t8_cmesh_new_periodic(comm, dim) - scheme = t8_scheme_new_default_cxx() + scheme = t8_scheme_new_default_cxx() - adapt_data = t8_step3_adapt_data_t( - (0.0, 0.0, 0.0), # Midpoints of the sphere. - 0.5, # Refine if inside this radius. - 0.7 # Coarsen if outside this radius. - ) + adapt_data = t8_step3_adapt_data_t((0.0, 0.0, 0.0), # Midpoints of the sphere. + 0.5, # Refine if inside this radius. + 0.7) - # Start with a uniform forest. - forest = t8_forest_new_uniform(cmesh, scheme, level, 0, comm) + # Start with a uniform forest. + forest = t8_forest_new_uniform(cmesh, scheme, level, 0, comm) - forest_apbg_ref = Ref(t8_forest_t()) - t8_forest_init(forest_apbg_ref) - forest_apbg = forest_apbg_ref[] + forest_apbg_ref = Ref(t8_forest_t()) + t8_forest_init(forest_apbg_ref) + forest_apbg = forest_apbg_ref[] - # Adapt, partition, balance and create ghost elements all in one go. - # See steps 3 and 4 for more details. - t8_forest_set_user_data(forest_apbg, Ref(adapt_data)) - t8_forest_set_adapt(forest_apbg, forest, @t8_adapt_callback(t8_step3_adapt_callback), 0) - t8_forest_set_partition(forest_apbg, C_NULL, 0) - t8_forest_set_balance(forest_apbg, C_NULL, 0) - t8_forest_set_ghost(forest_apbg, 1, T8_GHOST_FACES) - t8_forest_commit(forest_apbg) + # Adapt, partition, balance and create ghost elements all in one go. + # See steps 3 and 4 for more details. + t8_forest_set_user_data(forest_apbg, Ref(adapt_data)) + t8_forest_set_adapt(forest_apbg, forest, @t8_adapt_callback(t8_step3_adapt_callback), 0) + t8_forest_set_partition(forest_apbg, C_NULL, 0) + t8_forest_set_balance(forest_apbg, C_NULL, 0) + t8_forest_set_ghost(forest_apbg, 1, T8_GHOST_FACES) + t8_forest_commit(forest_apbg) - return forest_apbg + return forest_apbg end # Allocate and fill the element data array with `heights` from an arbitrary # mathematical function. Returns a pointer to the array which is then ownded by # the calling scope. function t8_step6_create_element_data(forest) - # Check that the forest is a committed. - @T8_ASSERT(t8_forest_is_committed(forest) == 1) + # Check that the forest is a committed. + @T8_ASSERT(t8_forest_is_committed(forest)==1) - # Get the number of local elements of forest. - num_local_elements = t8_forest_get_local_num_elements(forest) - # Get the number of ghost elements of forest. - num_ghost_elements = t8_forest_get_num_ghosts(forest) + # Get the number of local elements of forest. + num_local_elements = t8_forest_get_local_num_elements(forest) + # Get the number of ghost elements of forest. + num_ghost_elements = t8_forest_get_num_ghosts(forest) - # Build an array of our data that is as long as the number of elements plus - # the number of ghosts. - element_data = Array{t8_step6_data_per_element_t}(undef, num_local_elements + num_ghost_elements) + # Build an array of our data that is as long as the number of elements plus + # the number of ghosts. + element_data = Array{t8_step6_data_per_element_t}(undef, + num_local_elements + + num_ghost_elements) - # Get the number of trees that have elements of this process. - num_local_trees = t8_forest_get_num_local_trees(forest) + # Get the number of trees that have elements of this process. + num_local_trees = t8_forest_get_num_local_trees(forest) - # Compute vertex coordinates. Note: Julia has column-major. - verts = Matrix{Cdouble}(undef,3,3) + # Compute vertex coordinates. Note: Julia has column-major. + verts = Matrix{Cdouble}(undef, 3, 3) - # Element Midpoint - midpoint = Vector{Cdouble}(undef,3) + # Element Midpoint + midpoint = Vector{Cdouble}(undef, 3) - # Loop over all local trees in the forest. - current_index = 0 - for itree = 0:num_local_trees-1 - tree_class = t8_forest_get_tree_class(forest, itree) - eclass_scheme = t8_forest_get_eclass_scheme(forest, tree_class) + # Loop over all local trees in the forest. + current_index = 0 + for itree in 0:(num_local_trees - 1) + tree_class = t8_forest_get_tree_class(forest, itree) + eclass_scheme = t8_forest_get_eclass_scheme(forest, tree_class) - # Get the number of elements of this tree. - num_elements_in_tree = t8_forest_get_tree_num_elements(forest, itree) + # Get the number of elements of this tree. + num_elements_in_tree = t8_forest_get_tree_num_elements(forest, itree) - # Loop over all local elements in the tree. - for ielement = 0:num_elements_in_tree-1 - current_index += 1 # Note: Julia has 1-based indexing, while C/C++ starts with 0. + # Loop over all local elements in the tree. + for ielement in 0:(num_elements_in_tree - 1) + current_index += 1 # Note: Julia has 1-based indexing, while C/C++ starts with 0. - element = t8_forest_get_element_in_tree(forest, itree, ielement) + element = t8_forest_get_element_in_tree(forest, itree, ielement) - level = t8_element_level(eclass_scheme, element) - volume = t8_forest_element_volume(forest, itree, element) + level = t8_element_level(eclass_scheme, element) + volume = t8_forest_element_volume(forest, itree, element) - t8_forest_element_centroid(forest, itree, element, pointer(midpoint)) + t8_forest_element_centroid(forest, itree, element, pointer(midpoint)) - t8_element_vertex_reference_coords(eclass_scheme, element, 0, @view(verts[:,1])) - t8_element_vertex_reference_coords(eclass_scheme, element, 1, @view(verts[:,2])) - t8_element_vertex_reference_coords(eclass_scheme, element, 2, @view(verts[:,3])) + t8_element_vertex_reference_coords(eclass_scheme, element, 0, + @view(verts[:, 1])) + t8_element_vertex_reference_coords(eclass_scheme, element, 1, + @view(verts[:, 2])) + t8_element_vertex_reference_coords(eclass_scheme, element, 2, + @view(verts[:, 3])) - dx = verts[1,2] - verts[1,1] - dy = verts[2,3] - verts[2,1] + dx = verts[1, 2] - verts[1, 1] + dy = verts[2, 3] - verts[2, 1] - # Shift x and y to the center since the domain is [0,1] x [0,1]. - x = midpoint[1] - 0.5 - y = midpoint[2] - 0.5 - r = sqrt(x * x + y * y) * 20.0 # arbitrarily scaled radius + # Shift x and y to the center since the domain is [0,1] x [0,1]. + x = midpoint[1] - 0.5 + y = midpoint[2] - 0.5 + r = sqrt(x * x + y * y) * 20.0 # arbitrarily scaled radius - # Some 'interesting' height function. - height = sin(2.0 * r) / r + # Some 'interesting' height function. + height = sin(2.0 * r) / r - element_data[current_index] = t8_step6_data_per_element_t( - level, volume, Tuple(midpoint), dx, dy, height, 0.0, 0.0 - ) + element_data[current_index] = t8_step6_data_per_element_t(level, volume, + Tuple(midpoint), dx, + dy, height, 0.0, 0.0) + end end - end - return element_data + return element_data end # Gather the 3x3 stencil for each element and compute finite difference approximations # for schlieren and curvature of the stored heights in the elements. function t8_step6_compute_stencil(forest, element_data) - # Check that forest is a committed, that is valid and usable, forest. - @T8_ASSERT(t8_forest_is_committed(forest) == 1) - - # Get the number of trees that have elements of this process. - num_local_trees = t8_forest_get_num_local_trees(forest) - - stencil = Matrix{Cdouble}(undef, 3, 3) - dx = Vector{Cdouble}(undef, 3) - dy = Vector{Cdouble}(undef, 3) - - # Loop over all local trees in the forest. For each local tree the element - # data (level, midpoint[3], dx, dy, volume, height, schlieren, curvature) of - # each element is calculated and stored into the element data array. - current_index = 0 - for itree = 0:num_local_trees-1 - tree_class = t8_forest_get_tree_class(forest, itree) - eclass_scheme = t8_forest_get_eclass_scheme(forest, tree_class) - - num_elements_in_tree = t8_forest_get_tree_num_elements(forest, itree) - - # Loop over all local elements in the tree. - for ielement = 0:num_elements_in_tree-1 - current_index += 1 # Note: Julia has 1-based indexing, while C/C++ starts with 0. - - element = t8_forest_get_element_in_tree(forest, itree, ielement) - - # Gather center point of the 3x3 stencil. - stencil[2,2] = element_data[current_index].height - dx[2] = element_data[current_index].dx - dy[2] = element_data[current_index].dy - - # Loop over all faces of an element. - num_faces = t8_element_num_faces(eclass_scheme, element) - for iface = 1:num_faces - neighids_ref = Ref{Ptr{t8_locidx_t}}() - neighbors_ref = Ref{Ptr{Ptr{t8_element}}}() - neigh_scheme_ref = Ref{Ptr{t8_eclass_scheme}}() - - dual_faces_ref = Ref{Ptr{Cint}}() - num_neighbors_ref = Ref{Cint}() - - forest_is_balanced = Cint(1) - - t8_forest_leaf_face_neighbors(forest, itree, element, - neighbors_ref, iface-1, dual_faces_ref, num_neighbors_ref, - neighids_ref, neigh_scheme_ref, forest_is_balanced) - - num_neighbors = num_neighbors_ref[] - dual_faces = 1 .+ unsafe_wrap(Array, dual_faces_ref[], num_neighbors) - neighids = 1 .+ unsafe_wrap(Array, neighids_ref[], num_neighbors) - neighbors = unsafe_wrap(Array, neighbors_ref[], num_neighbors) - neigh_scheme = neigh_scheme_ref[] - - # Retrieve the `height` of the face neighbor. Account for two neighbors - # in case of a non-conforming interface by computing the average. - height = 0.0 - if num_neighbors > 0 - for ineigh = 1:num_neighbors - height = height + element_data[neighids[ineigh]].height - end - height = height / num_neighbors - end - - # Fill in the neighbor information of the 3x3 stencil. - if iface == 1 # NORTH - stencil[1,2] = height - dx[1] = element_data[neighids[1]].dx - elseif iface == 2 # SOUTH - stencil[3,2] = height - dx[3] = element_data[neighids[1]].dx - elseif iface == 3 # WEST - stencil[2,1] = height - dy[1] = element_data[neighids[1]].dy - elseif iface == 4 # EAST - stencil[2,3] = height - dy[3] = element_data[neighids[1]].dy + # Check that forest is a committed, that is valid and usable, forest. + @T8_ASSERT(t8_forest_is_committed(forest)==1) + + # Get the number of trees that have elements of this process. + num_local_trees = t8_forest_get_num_local_trees(forest) + + stencil = Matrix{Cdouble}(undef, 3, 3) + dx = Vector{Cdouble}(undef, 3) + dy = Vector{Cdouble}(undef, 3) + + # Loop over all local trees in the forest. For each local tree the element + # data (level, midpoint[3], dx, dy, volume, height, schlieren, curvature) of + # each element is calculated and stored into the element data array. + current_index = 0 + for itree in 0:(num_local_trees - 1) + tree_class = t8_forest_get_tree_class(forest, itree) + eclass_scheme = t8_forest_get_eclass_scheme(forest, tree_class) + + num_elements_in_tree = t8_forest_get_tree_num_elements(forest, itree) + + # Loop over all local elements in the tree. + for ielement in 0:(num_elements_in_tree - 1) + current_index += 1 # Note: Julia has 1-based indexing, while C/C++ starts with 0. + + element = t8_forest_get_element_in_tree(forest, itree, ielement) + + # Gather center point of the 3x3 stencil. + stencil[2, 2] = element_data[current_index].height + dx[2] = element_data[current_index].dx + dy[2] = element_data[current_index].dy + + # Loop over all faces of an element. + num_faces = t8_element_num_faces(eclass_scheme, element) + for iface in 1:num_faces + neighids_ref = Ref{Ptr{t8_locidx_t}}() + neighbors_ref = Ref{Ptr{Ptr{t8_element}}}() + neigh_scheme_ref = Ref{Ptr{t8_eclass_scheme}}() + + dual_faces_ref = Ref{Ptr{Cint}}() + num_neighbors_ref = Ref{Cint}() + + forest_is_balanced = Cint(1) + + t8_forest_leaf_face_neighbors(forest, itree, element, + neighbors_ref, iface - 1, dual_faces_ref, + num_neighbors_ref, + neighids_ref, neigh_scheme_ref, + forest_is_balanced) + + num_neighbors = num_neighbors_ref[] + dual_faces = 1 .+ unsafe_wrap(Array, dual_faces_ref[], num_neighbors) + neighids = 1 .+ unsafe_wrap(Array, neighids_ref[], num_neighbors) + neighbors = unsafe_wrap(Array, neighbors_ref[], num_neighbors) + neigh_scheme = neigh_scheme_ref[] + + # Retrieve the `height` of the face neighbor. Account for two neighbors + # in case of a non-conforming interface by computing the average. + height = 0.0 + if num_neighbors > 0 + for ineigh in 1:num_neighbors + height = height + element_data[neighids[ineigh]].height + end + height = height / num_neighbors + end + + # Fill in the neighbor information of the 3x3 stencil. + if iface == 1 # NORTH + stencil[1, 2] = height + dx[1] = element_data[neighids[1]].dx + elseif iface == 2 # SOUTH + stencil[3, 2] = height + dx[3] = element_data[neighids[1]].dx + elseif iface == 3 # WEST + stencil[2, 1] = height + dy[1] = element_data[neighids[1]].dy + elseif iface == 4 # EAST + stencil[2, 3] = height + dy[3] = element_data[neighids[1]].dy + end + + # Free allocated memory. + sc_free(t8_get_package_id(), neighbors_ref[]) + sc_free(t8_get_package_id(), dual_faces_ref[]) + sc_free(t8_get_package_id(), neighids_ref[]) + end + + # Prepare finite difference computations. The code also accounts for non-conforming interfaces. + xslope_m = 0.5 / (dx[1] + dx[2]) * (stencil[2, 2] - stencil[1, 2]) + xslope_p = 0.5 / (dx[2] + dx[3]) * (stencil[3, 2] - stencil[2, 2]) + + yslope_m = 0.5 / (dy[1] + dy[2]) * (stencil[2, 2] - stencil[2, 1]) + yslope_p = 0.5 / (dy[2] + dy[3]) * (stencil[2, 3] - stencil[2, 2]) + + xslope = 0.5 * (xslope_m + xslope_p) + yslope = 0.5 * (yslope_m + yslope_p) + + # TODO: Probably still not optimal at non-conforming interfaces. + xcurve = (xslope_p - xslope_m) * 4 / (dx[1] + 2.0 * dx[2] + dx[3]) + ycurve = (yslope_p - yslope_m) * 4 / (dy[1] + 2.0 * dy[2] + dy[3]) + + # Compute schlieren and curvature norm. + schlieren = sqrt(xslope * xslope + yslope * yslope) + curvature = sqrt(xcurve * xcurve + ycurve * ycurve) + + element_data[current_index] = t8_step6_data_per_element_t(element_data[current_index].level, + element_data[current_index].volume, + element_data[current_index].midpoint, + element_data[current_index].dx, + element_data[current_index].dy, + element_data[current_index].height, + schlieren, + curvature) end - - # Free allocated memory. - sc_free(t8_get_package_id(), neighbors_ref[]) - sc_free(t8_get_package_id(), dual_faces_ref[]) - sc_free(t8_get_package_id(), neighids_ref[]) - end - - # Prepare finite difference computations. The code also accounts for non-conforming interfaces. - xslope_m = 0.5 / (dx[1] + dx[2]) * (stencil[2,2] - stencil[1,2]) - xslope_p = 0.5 / (dx[2] + dx[3]) * (stencil[3,2] - stencil[2,2]) - - yslope_m = 0.5 / (dy[1] + dy[2]) * (stencil[2,2] - stencil[2,1]) - yslope_p = 0.5 / (dy[2] + dy[3]) * (stencil[2,3] - stencil[2,2]) - - xslope = 0.5 * (xslope_m + xslope_p) - yslope = 0.5 * (yslope_m + yslope_p) - - # TODO: Probably still not optimal at non-conforming interfaces. - xcurve = (xslope_p - xslope_m) * 4 / (dx[1] + 2.0 * dx[2] + dx[3]) - ycurve = (yslope_p - yslope_m) * 4 / (dy[1] + 2.0 * dy[2] + dy[3]) - - # Compute schlieren and curvature norm. - schlieren = sqrt(xslope * xslope + yslope * yslope) - curvature = sqrt(xcurve * xcurve + ycurve * ycurve) - - element_data[current_index] = t8_step6_data_per_element_t( - element_data[current_index].level, - element_data[current_index].volume, - element_data[current_index].midpoint, - element_data[current_index].dx, - element_data[current_index].dy, - element_data[current_index].height, - schlieren, - curvature - ) end - end end # Each process has computed the data entries for its local elements. In order @@ -299,15 +302,17 @@ end # entries of our element data array with the value on the process that owns the # corresponding element. */ function t8_step6_exchange_ghost_data(forest, element_data) - # t8_forest_ghost_exchange_data expects an sc_array (of length num_local_elements + num_ghosts). - # We wrap our data array to an sc_array. - sc_array_wrapper = sc_array_new_data(pointer(element_data), sizeof(t8_step6_data_per_element_t), length(element_data)) + # t8_forest_ghost_exchange_data expects an sc_array (of length num_local_elements + num_ghosts). + # We wrap our data array to an sc_array. + sc_array_wrapper = sc_array_new_data(pointer(element_data), + sizeof(t8_step6_data_per_element_t), + length(element_data)) - # Carry out the data exchange. The entries with indices > num_local_elements will get overwritten. - t8_forest_ghost_exchange_data(forest, sc_array_wrapper) + # Carry out the data exchange. The entries with indices > num_local_elements will get overwritten. + t8_forest_ghost_exchange_data(forest, sc_array_wrapper) - # Destroy the wrapper array. This will not free the data memory since we used sc_array_new_data. - sc_array_destroy(sc_array_wrapper) + # Destroy the wrapper array. This will not free the data memory since we used sc_array_new_data. + sc_array_destroy(sc_array_wrapper) end # Write the forest as vtu and also write the element's volumes in the file. @@ -318,52 +323,46 @@ end # We support two types: T8_VTK_SCALAR - One double per element. # and T8_VTK_VECTOR - Three doubles per element. function t8_step6_output_data_to_vtu(forest, element_data, prefix) - num_elements = length(element_data) - - # We need to allocate a new array to store the data on their own. - # These arrays have one entry per local element. - heights = Vector{Cdouble}(undef, num_elements) - schlieren = Vector{Cdouble}(undef, num_elements) - curvature = Vector{Cdouble}(undef, num_elements) - - # Copy the element's volumes from our data array to the output array. - for ielem = 1:num_elements - heights[ielem] = element_data[ielem].height - schlieren[ielem] = element_data[ielem].schlieren - curvature[ielem] = element_data[ielem].curvature - end - - # WARNING: This code hangs for Julia v1.8.* or older. Use at least Julia v1.9. - vtk_data = [ - t8_vtk_data_field_t( - T8_VTK_SCALAR, - NTuple{8192, Cchar}(rpad("height\0", 8192, ' ')), - pointer(heights), - ), - t8_vtk_data_field_t( - T8_VTK_SCALAR, - NTuple{8192, Cchar}(rpad("schlieren\0", 8192, ' ')), - pointer(schlieren), - ), - t8_vtk_data_field_t( - T8_VTK_SCALAR, - NTuple{8192, Cchar}(rpad("curvature\0", 8192, ' ')), - pointer(curvature), - ) - ] - - # The number of user defined data fields to write. - num_data = length(vtk_data) - - # Write user defined data to vtu file. - write_treeid = 1 - write_mpirank = 1 - write_level = 1 - write_element_id = 1 - write_ghosts = 0 - t8_forest_write_vtk_ext(forest, prefix, write_treeid, write_mpirank, - write_level, write_element_id, write_ghosts, - 0, 0, num_data, pointer(vtk_data)) + num_elements = length(element_data) + + # We need to allocate a new array to store the data on their own. + # These arrays have one entry per local element. + heights = Vector{Cdouble}(undef, num_elements) + schlieren = Vector{Cdouble}(undef, num_elements) + curvature = Vector{Cdouble}(undef, num_elements) + + # Copy the element's volumes from our data array to the output array. + for ielem in 1:num_elements + heights[ielem] = element_data[ielem].height + schlieren[ielem] = element_data[ielem].schlieren + curvature[ielem] = element_data[ielem].curvature + end + + # WARNING: This code hangs for Julia v1.8.* or older. Use at least Julia v1.9. + vtk_data = [ + t8_vtk_data_field_t(T8_VTK_SCALAR, + NTuple{8192, Cchar}(rpad("height\0", 8192, ' ')), + pointer(heights)), + t8_vtk_data_field_t(T8_VTK_SCALAR, + NTuple{8192, Cchar}(rpad("schlieren\0", 8192, ' ')), + pointer(schlieren)), + t8_vtk_data_field_t(T8_VTK_SCALAR, + NTuple{8192, Cchar}(rpad("curvature\0", 8192, ' ')), + pointer(curvature)), + ] + + # The number of user defined data fields to write. + num_data = length(vtk_data) + + # Write user defined data to vtu file. + write_treeid = 1 + write_mpirank = 1 + write_level = 1 + write_element_id = 1 + write_ghosts = 0 + t8_forest_write_vtk_ext(forest, prefix, write_treeid, write_mpirank, + write_level, write_element_id, write_ghosts, + 0, 0, num_data, pointer(vtk_data)) end # The prefix for our output files. diff --git a/examples/t8_tutorial_build_cmesh.jl b/examples/t8_tutorial_build_cmesh.jl index 7202606..25fad15 100644 --- a/examples/t8_tutorial_build_cmesh.jl +++ b/examples/t8_tutorial_build_cmesh.jl @@ -132,246 +132,246 @@ using T8code.Libt8: SC_LP_PRODUCTION # +---+---+ function t8_cmesh_new_periodic_hybrid_2d(comm) - # 1. Defining an array with all vertices. - # Just all vertices of all trees. partly duplicated. - vertices = [ - 0, 0, 0, # tree 0, triangle - 0.5, 0, 0, - 0.5, 0.5, 0, - 0, 0, 0, # tree 1, triangle - 0.5, 0.5, 0, - 0, 0.5, 0, - 0.5, 0, 0, # tree 2, quad - 1, 0, 0, - 0.5, 0.5, 0, - 1, 0.5, 0, - 0, 0.5, 0, # tree 3, quad - 0.5, 0.5, 0, - 0, 1, 0, - 0.5, 1, 0, - 0.5, 0.5, 0, # tree 4, triangle - 1, 0.5, 0, - 1, 1, 0, - 0.5, 0.5, 0, # tree 5, triangle - 1, 1, 0, - 0.5, 1, 0 - ] - - # 2. Initialization of the mesh. - cmesh_ref = Ref(t8_cmesh_t()) - t8_cmesh_init(cmesh_ref) - cmesh = cmesh_ref[] - - # 3. Definition of the geometry. - linear_geom = t8_geometry_linear_new(2) - t8_cmesh_register_geometry(cmesh, linear_geom) # Use linear geometry. - - # 4. Definition of the classes of the different trees. - t8_cmesh_set_tree_class(cmesh, 0, T8_ECLASS_TRIANGLE) - t8_cmesh_set_tree_class(cmesh, 1, T8_ECLASS_TRIANGLE) - t8_cmesh_set_tree_class(cmesh, 2, T8_ECLASS_QUAD) - t8_cmesh_set_tree_class(cmesh, 3, T8_ECLASS_QUAD) - t8_cmesh_set_tree_class(cmesh, 4, T8_ECLASS_TRIANGLE) - t8_cmesh_set_tree_class(cmesh, 5, T8_ECLASS_TRIANGLE) - - # 5. Classification of the vertices for each tree. - t8_cmesh_set_tree_vertices(cmesh, 0, vertices, 3) - t8_cmesh_set_tree_vertices(cmesh, 1, vertices + 9, 3) - t8_cmesh_set_tree_vertices(cmesh, 2, vertices + 18, 4) - t8_cmesh_set_tree_vertices(cmesh, 3, vertices + 30, 4) - t8_cmesh_set_tree_vertices(cmesh, 4, vertices + 42, 3) - t8_cmesh_set_tree_vertices(cmesh, 5, vertices + 51, 3) - - # 6. Definition of the face neighbors between the different trees. - t8_cmesh_set_join(cmesh, 0, 1, 1, 2, 0) - t8_cmesh_set_join(cmesh, 0, 2, 0, 0, 0) - t8_cmesh_set_join(cmesh, 0, 3, 2, 3, 0) - - t8_cmesh_set_join(cmesh, 1, 3, 0, 2, 1) - t8_cmesh_set_join(cmesh, 1, 2, 1, 1, 0) - - t8_cmesh_set_join(cmesh, 2, 4, 3, 2, 0) - t8_cmesh_set_join(cmesh, 2, 5, 2, 0, 1) - - t8_cmesh_set_join(cmesh, 3, 5, 1, 1, 0) - t8_cmesh_set_join(cmesh, 3, 4, 0, 0, 0) - - t8_cmesh_set_join(cmesh, 4, 5, 1, 2, 0) - - # 7. Commit the mesh. - t8_cmesh_commit(cmesh, comm) - - return cmesh + # 1. Defining an array with all vertices. + # Just all vertices of all trees. partly duplicated. + vertices = [ + 0, 0, 0, # tree 0, triangle + 0.5, 0, 0, + 0.5, 0.5, 0, + 0, 0, 0, # tree 1, triangle + 0.5, 0.5, 0, + 0, 0.5, 0, + 0.5, 0, 0, # tree 2, quad + 1, 0, 0, + 0.5, 0.5, 0, + 1, 0.5, 0, + 0, 0.5, 0, # tree 3, quad + 0.5, 0.5, 0, + 0, 1, 0, + 0.5, 1, 0, + 0.5, 0.5, 0, # tree 4, triangle + 1, 0.5, 0, + 1, 1, 0, + 0.5, 0.5, 0, # tree 5, triangle + 1, 1, 0, + 0.5, 1, 0, + ] + + # 2. Initialization of the mesh. + cmesh_ref = Ref(t8_cmesh_t()) + t8_cmesh_init(cmesh_ref) + cmesh = cmesh_ref[] + + # 3. Definition of the geometry. + linear_geom = t8_geometry_linear_new(2) + t8_cmesh_register_geometry(cmesh, linear_geom) # Use linear geometry. + + # 4. Definition of the classes of the different trees. + t8_cmesh_set_tree_class(cmesh, 0, T8_ECLASS_TRIANGLE) + t8_cmesh_set_tree_class(cmesh, 1, T8_ECLASS_TRIANGLE) + t8_cmesh_set_tree_class(cmesh, 2, T8_ECLASS_QUAD) + t8_cmesh_set_tree_class(cmesh, 3, T8_ECLASS_QUAD) + t8_cmesh_set_tree_class(cmesh, 4, T8_ECLASS_TRIANGLE) + t8_cmesh_set_tree_class(cmesh, 5, T8_ECLASS_TRIANGLE) + + # 5. Classification of the vertices for each tree. + t8_cmesh_set_tree_vertices(cmesh, 0, vertices, 3) + t8_cmesh_set_tree_vertices(cmesh, 1, vertices + 9, 3) + t8_cmesh_set_tree_vertices(cmesh, 2, vertices + 18, 4) + t8_cmesh_set_tree_vertices(cmesh, 3, vertices + 30, 4) + t8_cmesh_set_tree_vertices(cmesh, 4, vertices + 42, 3) + t8_cmesh_set_tree_vertices(cmesh, 5, vertices + 51, 3) + + # 6. Definition of the face neighbors between the different trees. + t8_cmesh_set_join(cmesh, 0, 1, 1, 2, 0) + t8_cmesh_set_join(cmesh, 0, 2, 0, 0, 0) + t8_cmesh_set_join(cmesh, 0, 3, 2, 3, 0) + + t8_cmesh_set_join(cmesh, 1, 3, 0, 2, 1) + t8_cmesh_set_join(cmesh, 1, 2, 1, 1, 0) + + t8_cmesh_set_join(cmesh, 2, 4, 3, 2, 0) + t8_cmesh_set_join(cmesh, 2, 5, 2, 0, 1) + + t8_cmesh_set_join(cmesh, 3, 5, 1, 1, 0) + t8_cmesh_set_join(cmesh, 3, 4, 0, 0, 0) + + t8_cmesh_set_join(cmesh, 4, 5, 1, 2, 0) + + # 7. Commit the mesh. + t8_cmesh_commit(cmesh, comm) + + return cmesh end # Definition of a three dimensional mesh with linear geometry. # The mesh consists of two tetrahedra, two prisms, one pyramid, and one hexahedron. function t8_cmesh_new_hybrid_gate_3d(comm) - vertices = Vector{Cdouble}(undef, 24) - linear_geom = t8_geometry_linear_new(3) - - # Initialization of the mesh. - cmesh_ref = Ref(t8_cmesh_t()) - t8_cmesh_init(cmesh_ref) - cmesh = cmesh_ref[] - - # Definition of the geometry. - t8_cmesh_register_geometry(cmesh, linear_geom) # Use linear geometry. - - # Defitition of the classes of the different trees. - t8_cmesh_set_tree_class(cmesh, 0, T8_ECLASS_TET) - t8_cmesh_set_tree_class(cmesh, 1, T8_ECLASS_TET) - t8_cmesh_set_tree_class(cmesh, 2, T8_ECLASS_PRISM) - t8_cmesh_set_tree_class(cmesh, 3, T8_ECLASS_PRISM) - t8_cmesh_set_tree_class(cmesh, 4, T8_ECLASS_PYRAMID) - t8_cmesh_set_tree_class(cmesh, 5, T8_ECLASS_HEX) - - # Classification of the vertices for each tree. - t8_cmesh_set_join(cmesh, 0, 2, 0, 4, 0) - t8_cmesh_set_join(cmesh, 1, 3, 0, 4, 0) - t8_cmesh_set_join(cmesh, 2, 5, 0, 0, 0) - t8_cmesh_set_join(cmesh, 3, 5, 1, 1, 0) - t8_cmesh_set_join(cmesh, 4, 5, 4, 2, 0) - - # - # Definition of the first tree - # - - # Tetrahedron 1 vertices. - vertices[1] = 0.43 - vertices[2] = 0 - vertices[3] = 2 - - vertices[4] = 0 - vertices[5] = 0 - vertices[6] = 1 - - vertices[7] = 0.86 - vertices[8] = -0.5 - vertices[9] = 1 - - vertices[10] = 0.86 - vertices[11] = 0.5 - vertices[12] = 1 - - # Classification of the vertices for the first tree. - t8_cmesh_set_tree_vertices(cmesh, 0, vertices, 4) - - # - # Definition of the second tree - # - - # Tetrahedron 2 vertices. - for itet = 0:2 - vertices[ 1 + itet] = vertices[ 1 + itet] + (itet == 0 ? 1 + 0.86 : 0) - vertices[ 4 + itet] = vertices[ 7 + itet] + (itet == 0 ? 1 : 0) - vertices[10 + itet] = vertices[10 + itet] + (itet == 0 ? 1 : 0) - end - - vertices[7] = 1 + 2 * 0.86 - vertices[8] = 0 - vertices[9] = 1 - - # Classification of the vertices for the second tree. - t8_cmesh_set_tree_vertices(cmesh, 1, vertices, 4) - - # - # Definition of the third tree - # - - # Prism 1 vertices. - vertices[1] = 0 - vertices[2] = 0 - vertices[3] = 0 - - vertices[4] = 0.86 - vertices[5] = -0.5 - vertices[6] = 0 - - vertices[7] = 0.86 - vertices[8] = 0.5 - vertices[9] = 0 - - # Translate +1 in z-axis for the upper vertices. - for iprism1 = 0:2 - vertices[10 + 3 * iprism1 ] = vertices[3 * iprism1 + 1] - vertices[10 + 3 * iprism1 + 1] = vertices[3 * iprism1 + 2] - vertices[10 + 3 * iprism1 + 2] = vertices[3 * iprism1 + 3] + 1 - end - - # Classification of the vertices for the third tree. - t8_cmesh_set_tree_vertices(cmesh, 2, vertices, 6) - - # - # Definition of the fourth tree - # - - # Prism 2 vertices. - for iprism2 = 0:2 - vertices[4 + iprism2] = vertices[1 + iprism2] + (iprism2 == 0 ? 1 + 2 * 0.86 : 0) - vertices[7 + iprism2] = vertices[7 + iprism2] + (iprism2 == 0 ? 1 : 0) - end - - vertices[1] = 0.86 + 1 - vertices[2] = -0.5 - vertices[3] = 0 - - # Translate +1 in z-axis for the upper vertices. - for iprism2 = 0:2 - vertices[10 + 3 * iprism2 ] = vertices[3 * iprism2 + 1] - vertices[10 + 3 * iprism2 + 1] = vertices[3 * iprism2 + 2] - vertices[10 + 3 * iprism2 + 2] = vertices[3 * iprism2 + 3] + 1 - end - - # Classification of the vertices for the fourth tree. - t8_cmesh_set_tree_vertices(cmesh, 3, vertices, 6) - - # - # Definition of the fifth tree - # - - # Pyramid vertices. - vertices[1] = 0.86 - vertices[2] = 0.5 - vertices[3] = 0 - - vertices[4] = 1.86 - vertices[5] = 0.5 - vertices[6] = 0 - - vertices[7] = 0.86 - vertices[8] = -0.5 - vertices[9] = 0 - - vertices[10] = 1.86 - vertices[11] = -0.5 - vertices[12] = 0 - - vertices[13] = 1.36 - vertices[14] = 0 - vertices[15] = -0.5 - - # Classification of the vertices for the fifth tree. - t8_cmesh_set_tree_vertices(cmesh, 4, vertices, 5) - - # - # Definition of the sixth tree - # - - # Hex coordinates. - for hex = 0:3 - vertices[3 * hex + 2] = vertices[3 * hex + 2] * (-1) - vertices[3 * hex + 13] = vertices[3 * hex + 1] - vertices[3 * hex + 14] = vertices[3 * hex + 2] - vertices[3 * hex + 15] = vertices[3 * hex + 3] + 1 - end - - # Classification of the vertices for the fifth tree. - t8_cmesh_set_tree_vertices(cmesh, 5, vertices, 8) - - # Commit the mesh. - t8_cmesh_commit(cmesh, comm) - return cmesh + vertices = Vector{Cdouble}(undef, 24) + linear_geom = t8_geometry_linear_new(3) + + # Initialization of the mesh. + cmesh_ref = Ref(t8_cmesh_t()) + t8_cmesh_init(cmesh_ref) + cmesh = cmesh_ref[] + + # Definition of the geometry. + t8_cmesh_register_geometry(cmesh, linear_geom) # Use linear geometry. + + # Defitition of the classes of the different trees. + t8_cmesh_set_tree_class(cmesh, 0, T8_ECLASS_TET) + t8_cmesh_set_tree_class(cmesh, 1, T8_ECLASS_TET) + t8_cmesh_set_tree_class(cmesh, 2, T8_ECLASS_PRISM) + t8_cmesh_set_tree_class(cmesh, 3, T8_ECLASS_PRISM) + t8_cmesh_set_tree_class(cmesh, 4, T8_ECLASS_PYRAMID) + t8_cmesh_set_tree_class(cmesh, 5, T8_ECLASS_HEX) + + # Classification of the vertices for each tree. + t8_cmesh_set_join(cmesh, 0, 2, 0, 4, 0) + t8_cmesh_set_join(cmesh, 1, 3, 0, 4, 0) + t8_cmesh_set_join(cmesh, 2, 5, 0, 0, 0) + t8_cmesh_set_join(cmesh, 3, 5, 1, 1, 0) + t8_cmesh_set_join(cmesh, 4, 5, 4, 2, 0) + + # + # Definition of the first tree + # + + # Tetrahedron 1 vertices. + vertices[1] = 0.43 + vertices[2] = 0 + vertices[3] = 2 + + vertices[4] = 0 + vertices[5] = 0 + vertices[6] = 1 + + vertices[7] = 0.86 + vertices[8] = -0.5 + vertices[9] = 1 + + vertices[10] = 0.86 + vertices[11] = 0.5 + vertices[12] = 1 + + # Classification of the vertices for the first tree. + t8_cmesh_set_tree_vertices(cmesh, 0, vertices, 4) + + # + # Definition of the second tree + # + + # Tetrahedron 2 vertices. + for itet in 0:2 + vertices[1 + itet] = vertices[1 + itet] + (itet == 0 ? 1 + 0.86 : 0) + vertices[4 + itet] = vertices[7 + itet] + (itet == 0 ? 1 : 0) + vertices[10 + itet] = vertices[10 + itet] + (itet == 0 ? 1 : 0) + end + + vertices[7] = 1 + 2 * 0.86 + vertices[8] = 0 + vertices[9] = 1 + + # Classification of the vertices for the second tree. + t8_cmesh_set_tree_vertices(cmesh, 1, vertices, 4) + + # + # Definition of the third tree + # + + # Prism 1 vertices. + vertices[1] = 0 + vertices[2] = 0 + vertices[3] = 0 + + vertices[4] = 0.86 + vertices[5] = -0.5 + vertices[6] = 0 + + vertices[7] = 0.86 + vertices[8] = 0.5 + vertices[9] = 0 + + # Translate +1 in z-axis for the upper vertices. + for iprism1 in 0:2 + vertices[10 + 3 * iprism1] = vertices[3 * iprism1 + 1] + vertices[10 + 3 * iprism1 + 1] = vertices[3 * iprism1 + 2] + vertices[10 + 3 * iprism1 + 2] = vertices[3 * iprism1 + 3] + 1 + end + + # Classification of the vertices for the third tree. + t8_cmesh_set_tree_vertices(cmesh, 2, vertices, 6) + + # + # Definition of the fourth tree + # + + # Prism 2 vertices. + for iprism2 in 0:2 + vertices[4 + iprism2] = vertices[1 + iprism2] + (iprism2 == 0 ? 1 + 2 * 0.86 : 0) + vertices[7 + iprism2] = vertices[7 + iprism2] + (iprism2 == 0 ? 1 : 0) + end + + vertices[1] = 0.86 + 1 + vertices[2] = -0.5 + vertices[3] = 0 + + # Translate +1 in z-axis for the upper vertices. + for iprism2 in 0:2 + vertices[10 + 3 * iprism2] = vertices[3 * iprism2 + 1] + vertices[10 + 3 * iprism2 + 1] = vertices[3 * iprism2 + 2] + vertices[10 + 3 * iprism2 + 2] = vertices[3 * iprism2 + 3] + 1 + end + + # Classification of the vertices for the fourth tree. + t8_cmesh_set_tree_vertices(cmesh, 3, vertices, 6) + + # + # Definition of the fifth tree + # + + # Pyramid vertices. + vertices[1] = 0.86 + vertices[2] = 0.5 + vertices[3] = 0 + + vertices[4] = 1.86 + vertices[5] = 0.5 + vertices[6] = 0 + + vertices[7] = 0.86 + vertices[8] = -0.5 + vertices[9] = 0 + + vertices[10] = 1.86 + vertices[11] = -0.5 + vertices[12] = 0 + + vertices[13] = 1.36 + vertices[14] = 0 + vertices[15] = -0.5 + + # Classification of the vertices for the fifth tree. + t8_cmesh_set_tree_vertices(cmesh, 4, vertices, 5) + + # + # Definition of the sixth tree + # + + # Hex coordinates. + for hex in 0:3 + vertices[3 * hex + 2] = vertices[3 * hex + 2] * (-1) + vertices[3 * hex + 13] = vertices[3 * hex + 1] + vertices[3 * hex + 14] = vertices[3 * hex + 2] + vertices[3 * hex + 15] = vertices[3 * hex + 3] + 1 + end + + # Classification of the vertices for the fifth tree. + t8_cmesh_set_tree_vertices(cmesh, 5, vertices, 8) + + # Commit the mesh. + t8_cmesh_commit(cmesh, comm) + return cmesh end # The prefix for our output files. diff --git a/src/T8code.jl b/src/T8code.jl index 527b4b6..8fcfca2 100644 --- a/src/T8code.jl +++ b/src/T8code.jl @@ -43,7 +43,6 @@ Returns the version of the underlying `t8code` library (*not* of T8code.jl). """ version() = VersionNumber(T8_VERSION_MAJOR, T8_VERSION_MINOR) - const T8CODE_UUID = UUID("d0cc0030-9a40-4274-8435-baadcfd54fa1") """ @@ -153,7 +152,8 @@ path_sc_library() = _PREFERENCE_LIBSC Returns `false` if a system-provided MPI installation is set via the MPIPreferences, but not a system-provided `t8code` installation. In this case, T8code.jl is not usable. """ -preferences_set_correctly() = !(_PREFERENCE_LIBT8 == "t8code_jll" && MPIPreferences.binary == "system") +preferences_set_correctly() = !(_PREFERENCE_LIBT8 == "t8code_jll" && + MPIPreferences.binary == "system") const T8_QUAD_MAXLEVEL = 30 const T8_HEX_MAXLEVEL = 19 @@ -165,11 +165,11 @@ const t8_hex_root_len = 1 << T8_HEX_MAXLEVEL @inline t8_hex_len(l) = 1 << (T8_HEX_MAXLEVEL - l) macro T8_ASSERT(q) - :( $(esc(q)) ? nothing : throw(AssertionError($(string(q)))) ) + :($(esc(q)) ? nothing : throw(AssertionError($(string(q))))) end function t8_free(ptr) - Libt8.sc_free(t8_get_package_id(), ptr) + Libt8.sc_free(t8_get_package_id(), ptr) end # typedef int (*t8_forest_adapt_t) (t8_forest_t forest, @@ -181,7 +181,9 @@ end # const int num_elements, # t8_element_t *elements[]); macro t8_adapt_callback(callback) - :( @cfunction($callback, Cint, (Ptr{t8_forest}, Ptr{t8_forest}, t8_locidx_t, t8_locidx_t, Ptr{t8_eclass_scheme}, Cint, Cint, Ptr{Ptr{t8_element}})) ) + :(@cfunction($callback, Cint, + (Ptr{t8_forest}, Ptr{t8_forest}, t8_locidx_t, t8_locidx_t, + Ptr{t8_eclass_scheme}, Cint, Cint, Ptr{Ptr{t8_element}}))) end # typedef void (*t8_forest_replace_t) (t8_forest_t forest_old, @@ -193,12 +195,14 @@ end # int num_incoming, # t8_locidx_t first_incoming); macro t8_replace_callback(callback) - :( @cfunction($callback, Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, t8_locidx_t, Ptr{Cvoid}, Cint, Cint, t8_locidx_t, Cint, t8_locidx_t)) ) + :(@cfunction($callback, Cvoid, + (Ptr{Cvoid}, Ptr{Cvoid}, t8_locidx_t, Ptr{Cvoid}, Cint, Cint, t8_locidx_t, + Cint, t8_locidx_t))) end function __init__() if !preferences_set_correctly() - @warn "System MPI version detected, but not a system t8code version. To make T8code.jl work, you need to set the preferences, see https://github.com/DLR-AMR/T8code.jl#using-a-custom-version-of-mpi-andor-t8code." + @warn "System MPI version detected, but not a system t8code version. To make T8code.jl work, you need to set the preferences, see https://github.com/DLR-AMR/T8code.jl#using-a-custom-version-of-mpi-andor-t8code." end end diff --git a/test/cmesh/test_readmshfile.jl b/test/cmesh/test_readmshfile.jl index 3653a7a..373f8de 100644 --- a/test/cmesh/test_readmshfile.jl +++ b/test/cmesh/test_readmshfile.jl @@ -1,138 +1,128 @@ function t8_supported_msh_file(cmesh) - # Description of the properties of the example msh-files. - number_elements = 4 - elem_type = T8_ECLASS_TRIANGLE - - vertex = [ - [0, 0], - [2, 0], - [4, 0], - [1, 2], - [3, 2], - [2, 4], - ] - - # 0-based indexing - elements = [ - [0, 1, 3], - [1, 4, 3], - [1, 2, 4], - [3, 4, 5], - ] - - face_neigh_elem = [ - [1, -1,-1], - [3, 0, 2], - [-1, 1, -1], - [-1, -1, 1], - ] - - @assert cmesh != C_NULL - - # Checks if the cmesh was committed. - @assert t8_cmesh_is_committed(cmesh) == 1 - - # `t8_cmesh_is_face_consistend` is not part of the public API. - # Checks for face consistency. - # @assert t8_cmesh_trees_is_face_consistend(cmesh, cmesh->trees) "Cmesh face consistency failed." - - # Checks if the number of elements was read correctly. - @test t8_cmesh_get_num_trees(cmesh) == number_elements - - # Number of local trees. - lnum_trees = t8_cmesh_get_num_local_trees(cmesh) - # Iterate through the local elements and check if they were read properly. - # All trees should be local to the master rank. - for ltree_it = 0:lnum_trees-1 - tree_class = t8_cmesh_get_tree_class(cmesh, ltree_it) - @test t8_eclass_compare(tree_class, elem_type) == 0 || "Element type in msh-file was read incorrectly." - - # Get pointer to the vertices of the tree. - vertices_ptr = t8_cmesh_get_tree_vertices(cmesh, ltree_it) - vertices = unsafe_wrap(Array,vertices_ptr,9) - # Checking the msh-files elements and nodes. - for i = 0:2 - # Checks if x and y coordinate of the nodes are not read correctly. - @test vertex[elements[ltree_it+1][i+1]+1][1] == vertices[3*i + 1] || - "x coordinate was read incorrectly." - @test vertex[elements[ltree_it+1][i+1]+1][2] == vertices[3*i + 2] || - "y coordinate was read incorrectly" - - # Checks whether the face neighbor elements are not read correctly. - ltree_id = t8_cmesh_get_face_neighbor(cmesh, ltree_it, i, C_NULL, C_NULL) - @test ltree_id == face_neigh_elem[ltree_it+1][i+1] || - "The face neighbor element in the example test file was not read correctly." + # Description of the properties of the example msh-files. + number_elements = 4 + elem_type = T8_ECLASS_TRIANGLE + + vertex = [ + [0, 0], + [2, 0], + [4, 0], + [1, 2], + [3, 2], + [2, 4], + ] + + # 0-based indexing + elements = [ + [0, 1, 3], + [1, 4, 3], + [1, 2, 4], + [3, 4, 5], + ] + + face_neigh_elem = [ + [1, -1, -1], + [3, 0, 2], + [-1, 1, -1], + [-1, -1, 1], + ] + + @assert cmesh != C_NULL + + # Checks if the cmesh was committed. + @assert t8_cmesh_is_committed(cmesh) == 1 + + # `t8_cmesh_is_face_consistend` is not part of the public API. + # Checks for face consistency. + # @assert t8_cmesh_trees_is_face_consistend(cmesh, cmesh->trees) "Cmesh face consistency failed." + + # Checks if the number of elements was read correctly. + @test t8_cmesh_get_num_trees(cmesh) == number_elements + + # Number of local trees. + lnum_trees = t8_cmesh_get_num_local_trees(cmesh) + # Iterate through the local elements and check if they were read properly. + # All trees should be local to the master rank. + for ltree_it in 0:(lnum_trees - 1) + tree_class = t8_cmesh_get_tree_class(cmesh, ltree_it) + @test t8_eclass_compare(tree_class, elem_type) == 0 || + "Element type in msh-file was read incorrectly." + + # Get pointer to the vertices of the tree. + vertices_ptr = t8_cmesh_get_tree_vertices(cmesh, ltree_it) + vertices = unsafe_wrap(Array, vertices_ptr, 9) + # Checking the msh-files elements and nodes. + for i in 0:2 + # Checks if x and y coordinate of the nodes are not read correctly. + @test vertex[elements[ltree_it + 1][i + 1] + 1][1] == vertices[3 * i + 1] || + "x coordinate was read incorrectly." + @test vertex[elements[ltree_it + 1][i + 1] + 1][2] == vertices[3 * i + 2] || + "y coordinate was read incorrectly" + + # Checks whether the face neighbor elements are not read correctly. + ltree_id = t8_cmesh_get_face_neighbor(cmesh, ltree_it, i, C_NULL, C_NULL) + @test ltree_id == face_neigh_elem[ltree_it + 1][i + 1] || + "The face neighbor element in the example test file was not read correctly." + end end - end - end @testset "readmshfile" begin + @testset "test_msh_file_vers2_ascii" begin + fileprefix = "cmesh/testfiles/test_msh_file_vers2_ascii" + filename = fileprefix * ".msh" - @testset "test_msh_file_vers2_ascii" begin - - fileprefix = "cmesh/testfiles/test_msh_file_vers2_ascii" - filename = fileprefix * ".msh" - - @assert isfile(filename) "File not found: " * filename - - cmesh = t8_cmesh_from_msh_file(fileprefix, 1, comm, 2, 0, 0) - - @assert cmesh != C_NULL "Could not read cmesh from ascii version 2, but should be able to: " * filename - - t8_supported_msh_file(cmesh) - t8_cmesh_destroy(Ref(cmesh)) - - end + @assert isfile(filename) "File not found: "*filename - @testset "test_msh_file_vers4_ascii" begin + cmesh = t8_cmesh_from_msh_file(fileprefix, 1, comm, 2, 0, 0) - fileprefix = "cmesh/testfiles/test_msh_file_vers4_ascii" - filename = fileprefix * ".msh" + @assert cmesh!=C_NULL "Could not read cmesh from ascii version 2, but should be able to: "*filename - @assert isfile(filename) "File not found: " * filename + t8_supported_msh_file(cmesh) + t8_cmesh_destroy(Ref(cmesh)) + end - cmesh = t8_cmesh_from_msh_file(fileprefix, 1, comm, 2, 0, 0) + @testset "test_msh_file_vers4_ascii" begin + fileprefix = "cmesh/testfiles/test_msh_file_vers4_ascii" + filename = fileprefix * ".msh" - @assert cmesh != C_NULL "Could not read cmesh from ascii version 4, but should be able to: " * filename + @assert isfile(filename) "File not found: "*filename - t8_cmesh_destroy(Ref(cmesh)) + cmesh = t8_cmesh_from_msh_file(fileprefix, 1, comm, 2, 0, 0) - end + @assert cmesh!=C_NULL "Could not read cmesh from ascii version 4, but should be able to: "*filename - @testset "test_msh_file_vers2_bin" begin + t8_cmesh_destroy(Ref(cmesh)) + end - fileprefix = "cmesh/testfiles/test_msh_file_vers2_bin" - filename = fileprefix * ".msh" + @testset "test_msh_file_vers2_bin" begin + fileprefix = "cmesh/testfiles/test_msh_file_vers2_bin" + filename = fileprefix * ".msh" - @assert isfile(filename) "File not found: " * filename + @assert isfile(filename) "File not found: "*filename - cmesh = t8_cmesh_from_msh_file(fileprefix, 1, comm, 2, 0, 0) + cmesh = t8_cmesh_from_msh_file(fileprefix, 1, comm, 2, 0, 0) - @assert cmesh == C_NULL "Expected fail of reading binary msh file v.2, but did not fail." + @assert cmesh==C_NULL "Expected fail of reading binary msh file v.2, but did not fail." - if cmesh != C_NULL - t8_cmesh_destroy(Ref(cmesh)) + if cmesh != C_NULL + t8_cmesh_destroy(Ref(cmesh)) + end end - end - - @testset "test_msh_file_vers4_bin" begin + @testset "test_msh_file_vers4_bin" begin + fileprefix = "cmesh/testfiles/test_msh_file_vers4_bin" + filename = fileprefix * ".msh" - fileprefix = "cmesh/testfiles/test_msh_file_vers4_bin" - filename = fileprefix * ".msh" + @assert isfile(filename) "File not found: "*filename - @assert isfile(filename) "File not found: " * filename + cmesh = t8_cmesh_from_msh_file(fileprefix, 1, comm, 2, 0, 0) - cmesh = t8_cmesh_from_msh_file(fileprefix, 1, comm, 2, 0, 0) + @assert cmesh==C_NULL "Expected fail of reading binary msh file v.4, but did not fail." - @assert cmesh == C_NULL "Expected fail of reading binary msh file v.4, but did not fail." - - if cmesh != C_NULL - t8_cmesh_destroy(Ref(cmesh)) + if cmesh != C_NULL + t8_cmesh_destroy(Ref(cmesh)) + end end - - end - end diff --git a/test/forest/test_element_volume.jl b/test/forest/test_element_volume.jl index 5d3da3b..287ee63 100644 --- a/test/forest/test_element_volume.jl +++ b/test/forest/test_element_volume.jl @@ -13,64 +13,63 @@ # elements and compute it element-specific. # \param[in] pyra A pyramid # \return The volume of the pyramid -function pyramid_control_volume(eclass_scheme :: Ptr{t8_eclass_scheme}, pyra :: Ptr{t8_element}) :: Cdouble - level = t8_element_level(eclass_scheme, pyra) - # Both pyramids and tets have 1/8th of the parents volume, if the shape does not switch. - control_volume = 1.0 / 3.0 / (1 << (level * 3)) +function pyramid_control_volume(eclass_scheme::Ptr{t8_eclass_scheme}, + pyra::Ptr{t8_element})::Cdouble + level = t8_element_level(eclass_scheme, pyra) + # Both pyramids and tets have 1/8th of the parents volume, if the shape does not switch. + control_volume = 1.0 / 3.0 / (1 << (level * 3)) - # Ancestors switch the shape. A tetrahedron has a 1/16th of its parents - # volume. For all levels we already divided the control-volume by 8, hence - # we divide it by 2 once. - # NOTE: T8code.jl does not support this operation. - # if (pyra->switch_shape_at_level > 0) - # control_volume /= 2 - # end + # Ancestors switch the shape. A tetrahedron has a 1/16th of its parents + # volume. For all levels we already divided the control-volume by 8, hence + # we divide it by 2 once. + # NOTE: T8code.jl does not support this operation. + # if (pyra->switch_shape_at_level > 0) + # control_volume /= 2 + # end - return control_volume + return control_volume end @testset "element volume" begin + epsilon = 1e-9 - epsilon = 1e-9 + for eclass in T8_ECLASS_ZERO:t8_eclass(T8_ECLASS_COUNT - 1) + for level in 0:4 + scheme = t8_scheme_new_default_cxx() + cmesh = t8_cmesh_new_hypercube(t8_eclass(eclass), comm, 0, 0, 0) + forest = t8_forest_new_uniform(cmesh, scheme, level, 0, comm) - for eclass = T8_ECLASS_ZERO:t8_eclass(T8_ECLASS_COUNT-1) + # Compute the global number of elements. + global_num_elements = t8_forest_get_global_num_elements(forest) - for level = 0:4 - scheme = t8_scheme_new_default_cxx() - cmesh = t8_cmesh_new_hypercube(t8_eclass(eclass), comm, 0, 0, 0) - forest = t8_forest_new_uniform(cmesh, scheme, level, 0, comm) + # Vertices have a volume of 0. + control_volume = (eclass == T8_ECLASS_VERTEX) ? 0.0 : + (1.0 / global_num_elements) - # Compute the global number of elements. - global_num_elements = t8_forest_get_global_num_elements(forest) + local_num_trees = t8_forest_get_num_local_trees(forest) - # Vertices have a volume of 0. - control_volume = (eclass == T8_ECLASS_VERTEX) ? 0.0 : (1.0 / global_num_elements) + # Iterate over all elements. + for itree in 0:(local_num_trees - 1) + tree_class = t8_forest_get_tree_class(forest, itree) + eclass_scheme = t8_forest_get_eclass_scheme(forest, tree_class) - local_num_trees = t8_forest_get_num_local_trees(forest) - - # Iterate over all elements. - for itree = 0:local_num_trees-1 - tree_class = t8_forest_get_tree_class(forest, itree) - eclass_scheme = t8_forest_get_eclass_scheme(forest, tree_class) - - tree_elements = t8_forest_get_tree_num_elements(forest, itree) - for ielement = 0:tree_elements-1 - element = t8_forest_get_element_in_tree(forest, itree, ielement) - volume = t8_forest_element_volume(forest, itree, element) - if eclass == T8_ECLASS_PYRAMID - shape_volume = pyramid_control_volume(eclass_scheme, element) - # WORKAROUND: See `function pyramid_control_volume` for a comment - # why the following test is split into two. - @test abs(shape_volume-volume) <= epsilon || abs(shape_volume/2-volume) <= epsilon - else - @test abs(volume-control_volume) <= epsilon - end - end # ielement - end # itree - - t8_forest_unref(Ref(forest)) - end # level - - end # eclass + tree_elements = t8_forest_get_tree_num_elements(forest, itree) + for ielement in 0:(tree_elements - 1) + element = t8_forest_get_element_in_tree(forest, itree, ielement) + volume = t8_forest_element_volume(forest, itree, element) + if eclass == T8_ECLASS_PYRAMID + shape_volume = pyramid_control_volume(eclass_scheme, element) + # WORKAROUND: See `function pyramid_control_volume` for a comment + # why the following test is split into two. + @test abs(shape_volume - volume) <= epsilon || + abs(shape_volume / 2 - volume) <= epsilon + else + @test abs(volume - control_volume) <= epsilon + end + end # ielement + end # itree + t8_forest_unref(Ref(forest)) + end # level + end # eclass end # testset diff --git a/test/runtests.jl b/test/runtests.jl index d1cf05c..f509e0e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,26 +6,26 @@ import MPIPreferences @info "Testing T8code.jl with" MPIPreferences.binary MPIPreferences.abi @time @testset "T8code.jl tests" begin - # For some weird reason, the MPI tests must come first since they fail - # otherwise with a custom MPI installation. - @time @testset "MPI" begin - # Do a dummy `@test true`: - # If the process errors out the testset would error out as well, - # cf. https://github.com/JuliaParallel/MPI.jl/pull/391 - @test true + # For some weird reason, the MPI tests must come first since they fail + # otherwise with a custom MPI installation. + @time @testset "MPI" begin + # Do a dummy `@test true`: + # If the process errors out the testset would error out as well, + # cf. https://github.com/JuliaParallel/MPI.jl/pull/391 + @test true - @info "Starting parallel tests" + @info "Starting parallel tests" - run(`$(mpiexec()) -n 2 $(Base.julia_cmd()) --threads=1 --check-bounds=yes --project=$(dirname(@__DIR__)) $(abspath("test_all.jl"))`) + run(`$(mpiexec()) -n 2 $(Base.julia_cmd()) --threads=1 --check-bounds=yes --project=$(dirname(@__DIR__)) $(abspath("test_all.jl"))`) - @info "Finished parallel tests" - end + @info "Finished parallel tests" + end - @time @testset "serial" begin - @info "Starting serial tests" + @time @testset "serial" begin + @info "Starting serial tests" - include("test_all.jl") + include("test_all.jl") - @info "Finished serial tests" - end + @info "Finished serial tests" + end end diff --git a/test/test_all.jl b/test/test_all.jl index 008b853..88f6f7e 100644 --- a/test/test_all.jl +++ b/test/test_all.jl @@ -13,22 +13,22 @@ MPI.Init() comm = MPI.COMM_WORLD @testset "init" begin - include("test_init.jl") + include("test_init.jl") end @testset "general" begin - if !Sys.iswindows() - # These tests do not work in Windows since the DLL loader does not search for symbols beyond libt8.dll. - include("test_refcount.jl") - end + if !Sys.iswindows() + # These tests do not work in Windows since the DLL loader does not search for symbols beyond libt8.dll. + include("test_refcount.jl") + end end @testset "cmesh" begin - include("cmesh/test_readmshfile.jl") + include("cmesh/test_readmshfile.jl") end @testset "forest" begin - include("forest/test_element_volume.jl") + include("forest/test_element_volume.jl") end end # module diff --git a/test/test_init.jl b/test/test_init.jl index 4e29ef0..46b2415 100644 --- a/test/test_init.jl +++ b/test/test_init.jl @@ -1,14 +1,14 @@ @testset "T8code.uses_mpi" begin - @test T8code.uses_mpi() == true + @test T8code.uses_mpi() == true end @testset "T8code.init" begin - @test_nowarn T8code.Libt8.sc_init(comm, 1, 1, C_NULL, SC_LP_DEFAULT) - @test_nowarn t8_init(SC_LP_DEFAULT) + @test_nowarn T8code.Libt8.sc_init(comm, 1, 1, C_NULL, SC_LP_DEFAULT) + @test_nowarn t8_init(SC_LP_DEFAULT) end @testset "sc_ functions" begin - @test_nowarn T8code.Libt8.sc_version() - @test_nowarn T8code.Libt8.sc_version_major() - @test_nowarn T8code.Libt8.sc_version_minor() + @test_nowarn T8code.Libt8.sc_version() + @test_nowarn T8code.Libt8.sc_version_major() + @test_nowarn T8code.Libt8.sc_version_minor() end diff --git a/test/test_refcount.jl b/test/test_refcount.jl index 8b3a3c7..067c1cd 100644 --- a/test/test_refcount.jl +++ b/test/test_refcount.jl @@ -1,79 +1,67 @@ @testset "test refcount" begin + @testset "test refcount init" begin + rc_ref = Ref(t8_refcount_t(-1, -1)) + t8_refcount_init(rc_ref) - @testset "test refcount init" begin + @test rc_ref[].refcount == 1 - rc_ref = Ref(t8_refcount_t(-1,-1)) - t8_refcount_init(rc_ref) - - @test rc_ref[].refcount == 1 - - @test T8code.Libt8.sc_refcount_is_active(rc_ref) == 1 - @test T8code.Libt8.sc_refcount_is_last(rc_ref) == 1 - @test T8code.Libt8.sc_refcount_unref(rc_ref) == 1 - - end - - @testset "test refcount init" begin - - rc_ptr = t8_refcount_new() - rc_ref = unsafe_wrap(Array,rc_ptr,1) - - @test rc_ref[].refcount == 1 - - @test T8code.Libt8.sc_refcount_is_active(rc_ptr) == 1 - @test T8code.Libt8.sc_refcount_is_last(rc_ptr) == 1 - @test T8code.Libt8.sc_refcount_unref(rc_ptr) == 1 - - t8_refcount_destroy(rc_ptr) - - end - - @testset "test refcount IsActive" begin - - rc_ref = Ref(t8_refcount_t(-1,-1)) - t8_refcount_init(rc_ref) - - @test rc_ref[].refcount == 1 + @test T8code.Libt8.sc_refcount_is_active(rc_ref) == 1 + @test T8code.Libt8.sc_refcount_is_last(rc_ref) == 1 + @test T8code.Libt8.sc_refcount_unref(rc_ref) == 1 + end - @test T8code.Libt8.sc_refcount_is_active(rc_ref) == 1 - @test T8code.Libt8.sc_refcount_unref(rc_ref) == 1 - @test T8code.Libt8.sc_refcount_is_active(rc_ref) == 0 + @testset "test refcount init" begin + rc_ptr = t8_refcount_new() + rc_ref = unsafe_wrap(Array, rc_ptr, 1) - end + @test rc_ref[].refcount == 1 - @testset "test refcount IsLast" begin + @test T8code.Libt8.sc_refcount_is_active(rc_ptr) == 1 + @test T8code.Libt8.sc_refcount_is_last(rc_ptr) == 1 + @test T8code.Libt8.sc_refcount_unref(rc_ptr) == 1 - rc_ref = Ref(t8_refcount_t(-1,-1)) - t8_refcount_init(rc_ref) + t8_refcount_destroy(rc_ptr) + end - @test T8code.Libt8.sc_refcount_is_last(rc_ref) == 1 - T8code.Libt8.sc_refcount_ref(rc_ref) - @test T8code.Libt8.sc_refcount_is_last(rc_ref) == 0 + @testset "test refcount IsActive" begin + rc_ref = Ref(t8_refcount_t(-1, -1)) + t8_refcount_init(rc_ref) - @test T8code.Libt8.sc_refcount_unref(rc_ref) == 0 - @test T8code.Libt8.sc_refcount_unref(rc_ref) == 1 + @test rc_ref[].refcount == 1 - end + @test T8code.Libt8.sc_refcount_is_active(rc_ref) == 1 + @test T8code.Libt8.sc_refcount_unref(rc_ref) == 1 + @test T8code.Libt8.sc_refcount_is_active(rc_ref) == 0 + end - @testset "test refcount RefUnref" begin + @testset "test refcount IsLast" begin + rc_ref = Ref(t8_refcount_t(-1, -1)) + t8_refcount_init(rc_ref) - rc_ref = Ref(t8_refcount_t(-1,-1)) - t8_refcount_init(rc_ref) + @test T8code.Libt8.sc_refcount_is_last(rc_ref) == 1 + T8code.Libt8.sc_refcount_ref(rc_ref) + @test T8code.Libt8.sc_refcount_is_last(rc_ref) == 0 - for value = 1:9 - @test rc_ref[].refcount == value - T8code.Libt8.sc_refcount_ref(rc_ref) + @test T8code.Libt8.sc_refcount_unref(rc_ref) == 0 + @test T8code.Libt8.sc_refcount_unref(rc_ref) == 1 end - @test rc_ref[].refcount == 10 + @testset "test refcount RefUnref" begin + rc_ref = Ref(t8_refcount_t(-1, -1)) + t8_refcount_init(rc_ref) - for value = 9:-1:1 - @test T8code.Libt8.sc_refcount_unref(rc_ref) == 0 - @test rc_ref[].refcount == value - end + for value in 1:9 + @test rc_ref[].refcount == value + T8code.Libt8.sc_refcount_ref(rc_ref) + end - @test T8code.Libt8.sc_refcount_unref(rc_ref) == 1 + @test rc_ref[].refcount == 10 - end + for value in 9:-1:1 + @test T8code.Libt8.sc_refcount_unref(rc_ref) == 0 + @test rc_ref[].refcount == value + end + @test T8code.Libt8.sc_refcount_unref(rc_ref) == 1 + end end