formatting and include desired qfunction

Added CMake dependencies on individual data files

CMakeLists.txt

@@ -870,11 +870,12 @@ add_dependencies(exec
 #  - https://cmake.org/Bug/view.php?id=8438
 
 # Add a target to copy the mfem data directory to the build directory
-add_custom_command(OUTPUT data_is_copied
-  COMMAND ${CMAKE_COMMAND} -E copy_directory ${PROJECT_SOURCE_DIR}/data data
-  COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
-  COMMENT "Copying the data directory ...")
-add_custom_target(copy_data DEPENDS data_is_copied)
+# Implementable as a single copy_directory_if_different command w/ CMake >= 3.26
+file(GLOB DATA_FILES CONFIGURE_DEPENDS ${PROJECT_SOURCE_DIR}/data/*)
+add_custom_target(copy_data
+  COMMAND ${CMAKE_COMMAND} -E make_directory data
+  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${DATA_FILES} data
+  COMMENT "Syncing the data directory ...")
 # Add 'copy_data' as a prerequisite for all executables, if the source and the
 # build directories are not the same.
 if (NOT ("${PROJECT_SOURCE_DIR}" STREQUAL "${PROJECT_BINARY_DIR}"))
@@ -1006,6 +1007,10 @@ install(FILES
 install(EXPORT ${PROJECT_NAME_UC}Targets
     DESTINATION ${INSTALL_CMAKE_DIR})
 
+# Install the data directory if present, i.e. if the copy_data target is built
+install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/data
+  DESTINATION ${MFEM_INSTALL_DIR} OPTIONAL)
+
 #-------------------------------------------------------------------------------
 # Create 'config.mk' from 'config.mk.in' for the build and install locations and
 # define install rules for 'config.mk' and 'test.mk'

fem/dfem/doperator.hpp

@@ -968,12 +968,17 @@ void DifferentiableOperator::AddIntegrator(
          // Quadrature point local derivative cache for each element, with data
          // layout:
          // [test_vdim, test_op_dim, trial_vdim, trial_op_dim, qp, num_entities].
-         derivative_qp_caches[derivative_id] = Vector(test_vdim * test_op_dim *
-                                                      trial_vdim *
-                                                      total_trial_op_dim * num_qp * num_entities);
+         size_t derivative_qp_size = test_vdim * test_op_dim *
+                                     trial_vdim *
+                                     total_trial_op_dim * num_qp * num_entities;
+         if (derivative_qp_caches.count(derivative_qp_size) == 0)
+         {
+            derivative_qp_caches[derivative_qp_size] = Vector(derivative_qp_size);
+         }
+
          // Create local references for MSVC lambda capture compatibility
          auto& fields_ref = this->fields;
-         auto& derivative_qp_caches_ref = this->derivative_qp_caches[derivative_id];
+         auto& derivative_qp_caches_ref = this->derivative_qp_caches[derivative_qp_size];
 
          // In each of the callbacks we're saving the derivatives in the quadrature point
          // caches. This trades memory with computational effort but also minimizes

miniapps/gslib/CMakeLists.txt

@@ -15,16 +15,14 @@ set(MESH_GF_FILES
     triple-pt-1.gf
     triple-pt-2.gf
 )
-
+# add target which keeps required mesh files in sync
 set(SRC_MESH_GF_FILES)
 foreach(MESH_GF_FILE ${MESH_GF_FILES})
   list(APPEND SRC_MESH_GF_FILES ${CMAKE_CURRENT_SOURCE_DIR}/${MESH_GF_FILE})
 endforeach()
-add_custom_command(OUTPUT data_is_copied
+add_custom_target(copy_miniapps_gslib_data
   COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_GF_FILES} .
-  COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
-  COMMENT "Copying gslib miniapps data files ...")
-add_custom_target(copy_miniapps_gslib_data DEPENDS data_is_copied)
+  COMMENT "Syncing gslib miniapps data files ...")
 
 if (MFEM_USE_GSLIB)
   add_mfem_miniapp(schwarz_ex1

miniapps/meshing/CMakeLists.txt

@@ -27,11 +27,9 @@ set(SRC_MESH_FILES)
 foreach(MESH_FILE ${MESH_FILES})
   list(APPEND SRC_MESH_FILES ${CMAKE_CURRENT_SOURCE_DIR}/${MESH_FILE})
 endforeach()
-add_custom_command(OUTPUT data_is_copied
+add_custom_target(copy_miniapps_meshing_data
   COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_FILES} .
-  COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
-  COMMENT "Copying meshing miniapps data files ...")
-add_custom_target(copy_miniapps_meshing_data DEPENDS data_is_copied)
+  COMMENT "Syncing meshing miniapps data files ...")
 
 add_mfem_miniapp(klein-bottle
   MAIN klein-bottle.cpp

miniapps/multidomain/CMakeLists.txt

@@ -19,11 +19,9 @@ set(SRC_MESH_FILES)
 foreach(MESH_FILE ${MESH_FILES})
   list(APPEND SRC_MESH_FILES ${CMAKE_CURRENT_SOURCE_DIR}/${MESH_FILE})
 endforeach()
-add_custom_command(OUTPUT data_is_copied
+add_custom_target(copy_miniapps_multidomain_data
   COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_FILES} .
-  COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
-  COMMENT "Copying multidomain miniapps data files ...")
-add_custom_target(copy_miniapps_multidomain_data DEPENDS data_is_copied)
+  COMMENT "Syncing multidomain miniapps data files ...")
 
 # Parallel apps.
 if (MFEM_USE_MPI)

miniapps/nurbs/CMakeLists.txt

@@ -9,33 +9,12 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
-set(MESH_FILES
-  cube-nurbs.mesh
-  ijk-hex-nurbs.mesh
-  plus-nurbs-permuted.mesh
-  plus-nurbs.mesh
-  square-nurbs.mesh
-  two-cubes-nurbs-autoedge.mesh
-  two-cubes-nurbs-rot.mesh
-  two-cubes-nurbs.mesh
-  two-squares-nurbs-autoedge.mesh
-  two-squares-nurbs-rot.mesh
-  two-squares-nurbs.mesh
-  3patch-nurbs.mesh
-  3patch-nurbs-flipedge.mesh
-)
-# Add a target to copy the mesh files from the source directory; used by sample
-# runs.
-set(SRC_MESH_FILES)
-foreach(MESH_FILE ${MESH_FILES})
-  list(APPEND SRC_MESH_FILES ${CMAKE_CURRENT_SOURCE_DIR}/meshes/${MESH_FILE})
-endforeach()
-add_custom_command(OUTPUT data_is_copied
+# add target which keeps required mesh files in sync
+file(GLOB SRC_MESH_FILES CONFIGURE_DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/meshes/*)
+add_custom_target(copy_miniapps_nurbs_data
   COMMAND ${CMAKE_COMMAND} -E make_directory meshes
-  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_FILES} meshes/
-  COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
-  COMMENT "Copying nurbs miniapps data files ...")
-add_custom_target(copy_miniapps_nurbs_data DEPENDS data_is_copied)
+  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_FILES} meshes
+  COMMENT "Syncing nurbs miniapps data directory ...")
 
 add_mfem_miniapp(nurbs_ex1
   MAIN nurbs_ex1.cpp

miniapps/tools/CMakeLists.txt

@@ -20,11 +20,9 @@ set(SRC_MESH_GF_FILES)
 foreach(MESH_FILE ${MESH_GF_FILES})
   list(APPEND SRC_MESH_GF_FILES ${CMAKE_CURRENT_SOURCE_DIR}/${MESH_FILE})
 endforeach()
-add_custom_command(OUTPUT data_is_copied
-  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_GF_FILES} ../gslib/.
-  COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
-  COMMENT "Copying tools miniapps data files ...")
-add_custom_target(copy_miniapps_tools_data DEPENDS data_is_copied)
+add_custom_target(copy_miniapps_tools_data
+  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_GF_FILES} .
+  COMMENT "Syncing tools miniapps data files ...")
 
 add_mfem_miniapp(display-basis
   MAIN display-basis.cpp

tests/unit/CMakeLists.txt

@@ -11,6 +11,24 @@
 
 project(mfem-unit-tests NONE)
 
+# Define a target that all examples and miniapps will depend on.
+set(MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME test_exec_prerequisites)
+add_custom_target(${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME})
+
+# Add a target to copy the mfem data directory to the build directory
+# Implementable as a single copy_directory_if_different command w/ CMake >= 3.26
+file(GLOB DATA_FILES CONFIGURE_DEPENDS ${PROJECT_SOURCE_DIR}/data/*)
+add_custom_target(copy_test_data
+  COMMAND ${CMAKE_COMMAND} -E make_directory data
+  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${DATA_FILES} data
+  COMMENT "Syncing the test data directory ...")
+# Add 'copy_test_data' as a prerequisite for test executables, if the source and the
+# build directories are not the same.
+if (NOT ("${PROJECT_SOURCE_DIR}" STREQUAL "${PROJECT_BINARY_DIR}"))
+  add_dependencies(${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME} copy_test_data)
+  add_dependencies(${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME} copy_data)
+endif()
+
 # Include the source directory for the unit tests - catch.hpp is there.
 include_directories(BEFORE ${CMAKE_CURRENT_SOURCE_DIR})
 
@@ -167,25 +185,19 @@ endif()
 
 add_library(unit_tests_srcs OBJECT ${UNIT_TESTS_SRCS})
 target_link_libraries(unit_tests_srcs PUBLIC mfem)
+add_dependencies(unit_tests_srcs
+      ${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME})
 
 # All serial non-device unit tests are built into a single executable,
 # 'unit_tests'.
 mfem_add_executable(unit_tests unit_test_main.cpp)
 target_link_libraries(unit_tests unit_tests_srcs)
 add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} unit_tests)
-# Unit tests need the ../../data directory.
-add_dependencies(unit_tests copy_data)
 # ParSubMesh tests need meshes in ../../miniapps/multidomain
 add_dependencies(unit_tests copy_miniapps_multidomain_data)
 # NURBS tests need meshes in ../../miniapps/nurbs
 add_dependencies(unit_tests copy_miniapps_nurbs_data)
 
-# Copy data to the build directory.
-add_custom_command(TARGET unit_tests POST_BUILD
-COMMAND ${CMAKE_COMMAND} -E copy_directory
-        ${CMAKE_CURRENT_SOURCE_DIR}/data data
-        COMMENT "Copying the unit tests data directory ...")
-
 # Create a test called 'unit_tests' that runs the 'unit_tests' executable.
 # The unit tests can be built and run separately from the rest of the tests:
 #   make unit_tests
@@ -204,7 +216,6 @@ if (MFEM_USE_CUDA)
     set_property(SOURCE ${GPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
     mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS})
     target_link_libraries(gpu_unit_tests unit_tests_srcs)
-    add_dependencies(gpu_unit_tests copy_data)
     add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
     if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
         add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
@@ -220,7 +231,6 @@ if (MFEM_USE_HIP)
     set(GPU_UNIT_TESTS_SRCS gpu_unit_test_main.cpp)
     mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS})
     target_link_libraries(gpu_unit_tests unit_tests_srcs)
-    add_dependencies(gpu_unit_tests copy_data)
     add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
     if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
         add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
@@ -247,6 +257,8 @@ function(add_serial_miniapp_test name test_uvm)
     endif(MFEM_USE_HIP)
 
     mfem_add_executable(${name}_tests_cpu ${${NAME}_TESTS_SRCS})
+    add_dependencies(${name}_tests_cpu copy_miniapps_meshing_data)
+    add_dependencies(${name}_tests_cpu ${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME})
     target_compile_definitions(${name}_tests_cpu PUBLIC MFEM_${NAME}_DEVICE="cpu")
     target_link_libraries(${name}_tests_cpu mfem)
     add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ${name}_tests_cpu)
@@ -255,6 +267,8 @@ function(add_serial_miniapp_test name test_uvm)
     endif()
 
     mfem_add_executable(${name}_tests_debug ${${NAME}_TESTS_SRCS})
+    add_dependencies(${name}_tests_debug copy_miniapps_meshing_data)
+    add_dependencies(${name}_tests_debug ${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME})
     target_compile_definitions(${name}_tests_debug PUBLIC MFEM_${NAME}_DEVICE="debug")
     target_link_libraries(${name}_tests_debug mfem)
     add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ${name}_tests_debug)
@@ -264,6 +278,8 @@ function(add_serial_miniapp_test name test_uvm)
 
     if (MFEM_USE_CUDA OR MFEM_USE_HIP)
         mfem_add_executable(${name}_tests_gpu ${${NAME}_TESTS_SRCS})
+        add_dependencies(${name}_tests_gpu copy_miniapps_meshing_data)
+        add_dependencies(${name}_tests_gpu ${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME})
         target_compile_definitions(${name}_tests_gpu PUBLIC MFEM_${NAME}_DEVICE="gpu")
         target_link_libraries(${name}_tests_gpu mfem)
         add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ${name}_tests_gpu)
@@ -273,6 +289,8 @@ function(add_serial_miniapp_test name test_uvm)
 
         if (test_uvm)
             mfem_add_executable(${name}_tests_gpu_uvm ${${NAME}_TESTS_SRCS})
+            add_dependencies(${name}_tests_gpu_uvm copy_miniapps_meshing_data)
+            add_dependencies(${name}_tests_gpu_uvm ${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME})
             target_compile_definitions(${name}_tests_gpu_uvm PUBLIC
                MFEM_${NAME}_DEVICE="gpu:uvm")
             target_link_libraries(${name}_tests_gpu_uvm mfem)
@@ -288,7 +306,6 @@ endfunction(add_serial_miniapp_test)
 add_serial_miniapp_test(sedov ON)  # UVM ON
 add_serial_miniapp_test(tmop_pa OFF)  # UVM OFF
 # TMOP tests need meshes in ../../miniapps/meshing
-add_dependencies(tmop_pa_tests_cpu copy_miniapps_meshing_data)
 
 #-----------------------------------------------------------
 # SERIAL CEED TESTS:
@@ -342,7 +359,6 @@ if (MFEM_USE_MPI)
         set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
         set_property(SOURCE ${PGPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
         mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS})
-        add_dependencies(pgpu_unit_tests copy_data)
         target_link_libraries(pgpu_unit_tests unit_tests_srcs)
         add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
         foreach(np 1 ${MFEM_MPI_NP})
@@ -358,7 +374,6 @@ if (MFEM_USE_MPI)
         # pgpu_unit_tests
         set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
         mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS})
-        add_dependencies(pgpu_unit_tests copy_data)
         target_link_libraries(pgpu_unit_tests unit_tests_srcs)
         add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
         foreach(np 1 ${MFEM_MPI_NP})

tests/unit/dfem/test_diffusion.cpp

@@ -69,9 +69,9 @@ template <int DIM> struct Diffusion
 };
 
 template <int DIM>
-void diffusion(const char *filename, int p)
+void diffusion(const char *filename, int p, bool include_mass)
 {
-   CAPTURE(filename, DIM, p);
+   CAPTURE(filename, DIM, p, include_mass);
 
    Mesh smesh(filename);
    ParMesh pmesh(MPI_COMM_WORLD, smesh);
@@ -116,6 +116,10 @@ void diffusion(const char *filename, int p)
 
    ParBilinearForm blf_fa(&pfes);
    blf_fa.AddDomainIntegrator(new DiffusionIntegrator(rho_coeff, ir));
+   if(include_mass)
+   {
+      blf_fa.AddDomainIntegrator(new MassIntegrator(ir));
+   }
    blf_fa.SetAssemblyLevel(AssemblyLevel::FULL);
    blf_fa.Assemble();
    blf_fa.Finalize();
@@ -130,6 +134,16 @@ void diffusion(const char *filename, int p)
    static constexpr int U = 0, Coords = 1, Rho = 3;
    const auto sol = std::vector{ FieldDescriptor{ U, &pfes } };
 
+   // Define mass qfunction in case it's used.
+      const auto mf_mass_qf =
+         [] MFEM_HOST_DEVICE(const dscalar_t u, const real_t &rho /*Needed for some reason*/,
+                             const tensor<real_t, DIM, DIM> &J, const real_t &w)
+      { return tuple{u * w * det(J)}; };
+      const auto desired_mass_qf =
+         [] MFEM_HOST_DEVICE(const dscalar_t u,
+                             const tensor<real_t, DIM, DIM> &J, const real_t &w)
+      { return tuple{u * w * det(J)}; };
+
    SECTION("action")
    {
       DOperator dop_mf(sol, {{Rho, &rho_ps}, {Coords, mfes}}, pmesh);
@@ -139,6 +153,21 @@ void diffusion(const char *filename, int p)
                                         Gradient<Coords>{}, Weight{} },
                                  tuple{ Gradient<U>{} }, *ir,
                                  all_domain_attr);
+      if(include_mass)
+      {
+         // This version works, but had to include unused Rho parameter in qfunction.
+         dop_mf.AddDomainIntegrator(mf_mass_qf,
+                                    tuple{ Value<U>{}, Identity<Rho>{},
+                                           Gradient<Coords>{}, Weight{} },
+                                    tuple{ Value<U>{} }, *ir,
+                                    all_domain_attr);
+         // Below version does not work.
+         // dop_mf.AddDomainIntegrator(desired_mass_qf,
+         //                            tuple{ Value<U>{},
+         //                                   Gradient<Coords>{}, Weight{} },
+         //                            tuple{ Value<U>{} }, *ir,
+         //                            all_domain_attr);
+      }
       dop_mf.SetParameters({ &rho_coeff_cv, nodes });
 
       pfes.GetRestrictionMatrix()->Mult(x, X);
@@ -159,6 +188,9 @@ void diffusion(const char *filename, int p)
 
    SECTION("action partial assembly")
    {
+      // is there a better way to skip a test than this?
+      if(!include_mass)
+      {
       static constexpr int QData = 2;
       UniformParameterSpace qd_ps(pmesh, *ir, DIM * DIM);
       ParameterFunction qdata(qd_ps);
@@ -196,6 +228,7 @@ void diffusion(const char *filename, int p)
 
       REQUIRE(norm_global == MFEM_Approx(0.0));
       MPI_Barrier(MPI_COMM_WORLD);
+      }
    }
 
    SECTION("action linearized")
@@ -208,6 +241,14 @@ void diffusion(const char *filename, int p)
                                         Gradient<Coords>{}, Weight{} },
                                  tuple{ Gradient<U>{} }, *ir,
                                  all_domain_attr, derivatives);
+      if(include_mass)
+      {
+         dop_mf.AddDomainIntegrator(mf_mass_qf,
+                                    tuple{ Value<U>{}, Identity<Rho>{} /*need this along with dummy argument to q function for some reason*/,
+                                           Gradient<Coords>{}, Weight{} },
+                                    tuple{ Value<U>{} }, *ir,
+                                    all_domain_attr, derivatives);
+      }
       dop_mf.SetParameters({ &rho_coeff_cv, nodes });
       auto dRdU = dop_mf.GetDerivative(U, {&x}, {&rho_coeff_cv, nodes});
 
@@ -258,6 +299,10 @@ void diffusion(const char *filename, int p)
          ConstantCoefficient one(1.0);
          ParBilinearForm vblf_fa(&vpfes);
          vblf_fa.AddDomainIntegrator(new VectorDiffusionIntegrator(one, ir));
+         // if(include_mass)
+         // {
+         //    vblf_fa.AddDomainIntegrator(new VectorMassIntegrator(one, ir));
+         // }
          vblf_fa.SetAssemblyLevel(AssemblyLevel::LEGACYFULL);
          vblf_fa.Assemble();
          vblf_fa.Finalize();
@@ -290,6 +335,14 @@ void diffusion(const char *filename, int p)
                                         Gradient<Coords>{}, Weight{} },
                                  tuple{ Gradient<U>{} }, *ir,
                                  all_domain_attr, derivatives);
+      if(include_mass)
+      {
+         dop_mf.AddDomainIntegrator(mf_mass_qf,
+                                    tuple{ Value<U>{}, Identity<Rho>{},
+                                           Gradient<Coords>{}, Weight{} },
+                                    tuple{ Value<U>{} }, *ir,
+                                    all_domain_attr, derivatives);
+      }
       dop_mf.SetParameters({ &rho_coeff_cv, nodes });
       auto dRdU = dop_mf.GetDerivative(U, {&x}, {&rho_coeff_cv, nodes});
 
@@ -306,6 +359,8 @@ TEST_CASE("dFEM Diffusion", "[Parallel][dFEM][GPU]")
 
    const auto p = !all_tests ? 2 : GENERATE(1, 2, 3);
 
+   const bool include_mass = GENERATE(false, true);
+
    SECTION("2d")
    {
       const auto filename2d =
@@ -316,7 +371,7 @@ TEST_CASE("dFEM Diffusion", "[Parallel][dFEM][GPU]")
             "../../data/inline-quad.mesh",
             "../../data/periodic-square.mesh"
          );
-      diffusion<2>(filename2d, p);
+      diffusion<2>(filename2d, p, include_mass);
    }
 
    SECTION("3d")
@@ -329,7 +384,7 @@ TEST_CASE("dFEM Diffusion", "[Parallel][dFEM][GPU]")
             "../../data/toroid-hex.mesh",
             "../../data/periodic-cube.mesh"
          );
-      diffusion<3>(filename3d, p);
+      diffusion<3>(filename3d, p, include_mass);
    }
 }