diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index aab1e9e636889af11d003a79446fa064450e5be1..fc71924b38c7740bae5f0df4fb53f23dec048fd7 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -50,17 +50,18 @@ void readGeometry(std::string const& fname, GeoLib::GEOObjects & geo_objects)
ProjectData::ProjectData() = default;
ProjectData::ProjectData(BaseLib::ConfigTree const& project_config,
- std::string const& path)
+ std::string const& project_directory,
+ std::string const& output_directory)
{
// geometry
std::string const geometry_file = BaseLib::copyPathToFileName(
- project_config.getConfParam<std::string>("geometry"), path
+ project_config.getConfParam<std::string>("geometry"), project_directory
);
detail::readGeometry(geometry_file, *_geoObjects);
// mesh
std::string const mesh_file = BaseLib::copyPathToFileName(
- project_config.getConfParam<std::string>("mesh"), path
+ project_config.getConfParam<std::string>("mesh"), project_directory
);
MeshLib::Mesh* const mesh = FileIO::readMeshFromFile(mesh_file);
@@ -82,7 +83,7 @@ ProjectData::ProjectData(BaseLib::ConfigTree const& project_config,
parseProcesses(project_config.getConfSubtree("processes"));
// output
- parseOutput(project_config.getConfSubtree("output"), path);
+ parseOutput(project_config.getConfSubtree("output"), output_directory);
// timestepping
parseTimeStepping(project_config.getConfSubtree("time_stepping"));
@@ -296,14 +297,12 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config)
}
void ProjectData::parseOutput(BaseLib::ConfigTree const& output_config,
- std::string const& path)
+ std::string const& output_directory)
{
output_config.checkConfParam("type", "VTK");
DBUG("Parse output configuration:");
- auto const file = output_config.getConfParam<std::string>("file");
-
- _output_file_prefix = path + file;
+ _output = ProcessLib::Output<GlobalSetupType>::newInstance(output_config, output_directory);
}
void ProjectData::parseTimeStepping(BaseLib::ConfigTree const& timestepping_config)
diff --git a/Applications/ApplicationsLib/ProjectData.h b/Applications/ApplicationsLib/ProjectData.h
index 87ca99dfc66e45eb34bcf1531e3a871b720e6ed8..2082ed79d5f2b680946d1b05e44ed07a6385c496 100644
--- a/Applications/ApplicationsLib/ProjectData.h
+++ b/Applications/ApplicationsLib/ProjectData.h
@@ -22,6 +22,7 @@
#include "ProcessLib/ProcessVariable.h"
#include "ProcessLib/Process.h"
#include "ProcessLib/Parameter.h"
+#include "ProcessLib/Output.h"
#include "NumLib/ODESolver/Types.h"
@@ -60,10 +61,13 @@ public:
ProjectData();
/// Constructs project data by parsing provided configuration.
- /// The additional path is used to find files referenced in the
- /// configuration.
+ ///
+ /// \param config_tree Configuration as read from the prj file.
+ /// \param project_directory Where to look for files referenced in the \c config_tree.
+ /// \param output_directory Where to write simulation output files to.
ProjectData(BaseLib::ConfigTree const& config_tree,
- std::string const& path);
+ std::string const& project_directory,
+ std::string const& output_directory);
ProjectData(ProjectData&) = delete;
@@ -131,10 +135,16 @@ public:
return _processes.end();
}
- std::string const&
- getOutputFilePrefix() const
+ ProcessLib::Output<GlobalSetupType> const&
+ getOutputControl() const
{
- return _output_file_prefix;
+ return *_output;
+ }
+
+ ProcessLib::Output<GlobalSetupType>&
+ getOutputControl()
+ {
+ return *_output;
}
TimeLoop& getTimeLoop()
@@ -174,7 +184,8 @@ private:
/// Parses the output configuration.
/// Parses the file tag and sets output file prefix.
- void parseOutput(BaseLib::ConfigTree const& output_config, std::string const& path);
+ void parseOutput(BaseLib::ConfigTree const& output_config,
+ std::string const& output_directory);
void parseTimeStepping(BaseLib::ConfigTree const& timestepping_config);
@@ -195,8 +206,7 @@ private:
/// Buffer for each parameter config passed to the process.
std::vector<std::unique_ptr<ProcessLib::ParameterBase>> _parameters;
- /// Output file path with project prefix.
- std::string _output_file_prefix;
+ std::unique_ptr<ProcessLib::Output<GlobalSetupType>> _output;
/// The time loop used to solve this project's processes.
std::unique_ptr<TimeLoop> _time_loop;
diff --git a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h
index 2223b22953470a9da1203df333f7b02fdddd373e..d2cccfbcc8d6aaa8b9d7a7205c69f0714e024038 100644
--- a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h
+++ b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h
@@ -29,7 +29,7 @@ template<typename Matrix, typename Vector>
class UncoupledProcessesTimeLoop
{
public:
- bool loop(ProjectData& project, std::string const& outdir);
+ bool loop(ProjectData& project);
~UncoupledProcessesTimeLoop();
@@ -151,11 +151,9 @@ private:
//! Solves one timestep for the given \c process.
bool solveOneTimeStepOneProcess(
- unsigned const timestep,
Vector& x, double const t, double const delta_t,
SingleProcessData& process_data,
- Process& process,
- std::string const& output_file_name);
+ Process& process);
//! Sets the EquationSystem for the given nonlinear solver,
//! which is Picard or Newton depending on the NLTag.
@@ -290,11 +288,9 @@ template<typename Matrix, typename Vector>
bool
UncoupledProcessesTimeLoop<Matrix, Vector>::
solveOneTimeStepOneProcess(
- unsigned const timestep,
Vector& x, double const t, double const delta_t,
SingleProcessData& process_data,
- typename UncoupledProcessesTimeLoop<Matrix, Vector>::Process& process,
- std::string const& output_file_name)
+ typename UncoupledProcessesTimeLoop<Matrix, Vector>::Process& process)
{
auto& time_disc = process.getTimeDiscretization();
auto& ode_sys = *process_data.tdisc_ode_sys;
@@ -314,7 +310,6 @@ solveOneTimeStepOneProcess(
time_disc.pushState(t, x, mat_strg);
process.postTimestep(x);
- process.output(output_file_name, timestep, x);
return nonlinear_solver_succeeded;
}
@@ -323,54 +318,79 @@ solveOneTimeStepOneProcess(
template<typename Matrix, typename Vector>
bool
UncoupledProcessesTimeLoop<Matrix, Vector>::
-loop(ProjectData& project, std::string const& outdir)
+loop(ProjectData& project)
{
auto per_process_data = initInternalData(project);
+ auto& out_ctrl = project.getOutputControl();
+ out_ctrl.init(project.processesBegin(), project.processesEnd());
+
auto const t0 = 0.0; // time of the IC
// init solution storage
setInitialConditions(project, t0, per_process_data);
+ // output initial conditions
+ {
+ unsigned pcs_idx = 0;
+ for (auto p = project.processesBegin(); p != project.processesEnd();
+ ++p, ++pcs_idx)
+ {
+ auto const& x0 = *_process_solutions[pcs_idx];
+ out_ctrl.doOutput(**p, 0, t0, x0);
+ }
+ }
+
auto const delta_t = 1.0;
// TODO only for now
// Make sure there will be exactly one iteration of the loop below.
- auto const t_end = 1.5;
+ auto const t_end = 0.5;
- double t;
+ double t, t_last;
unsigned timestep = 1; // the first timestep really is number one
bool nonlinear_solver_succeeded = true;
- for (t=t0+delta_t; t<t_end+delta_t; t+=delta_t, ++timestep)
+
+ for (t=t0+delta_t, t_last=t0;
+ t<t_end+delta_t;
+ t_last=t, t+=delta_t, ++timestep)
{
+ // TODO use process name
unsigned pcs_idx = 0;
for (auto p = project.processesBegin(); p != project.processesEnd();
++p, ++pcs_idx)
{
- std::string const& outpref = project.getOutputFilePrefix();
- std::string const output_file_name =
- BaseLib::joinPaths(outdir, outpref)
- + "_pcs_" + std::to_string(pcs_idx)
- + "_ts_" + std::to_string(timestep)
- // + "_t_" + std::to_string(t) // TODO: add that later
- + ".vtu";
-
auto& x = *_process_solutions[pcs_idx];
nonlinear_solver_succeeded = solveOneTimeStepOneProcess(
- timestep, x, t, delta_t,
- per_process_data[pcs_idx], **p, output_file_name);
+ x, t, delta_t, per_process_data[pcs_idx], **p);
if (!nonlinear_solver_succeeded) {
ERR("The nonlinear solver failed in timestep #%u at t = %g s"
" for process #%u.", timestep, t, pcs_idx);
+
+ // save unsuccessful solution
+ out_ctrl.doOutputAlways(**p, timestep, t, x);
+
break;
+ } else {
+ out_ctrl.doOutput(**p, timestep, t, x);
}
}
if (!nonlinear_solver_succeeded) break;
+ }
- break; // TODO only do a single timestep for now
+ // output last timestep
+ if (nonlinear_solver_succeeded)
+ {
+ unsigned pcs_idx = 0;
+ for (auto p = project.processesBegin(); p != project.processesEnd();
+ ++p, ++pcs_idx)
+ {
+ auto const& x = *_process_solutions[pcs_idx];
+ out_ctrl.doOutputLastTimestep(**p, timestep-1, t_last, x);
+ }
}
return nonlinear_solver_succeeded;
diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake
index a9f658c5465b2bbd0d9f968da26fefb3406e10d5..49533168851b90000576b2d19fbc0ee53cd46588 100644
--- a/Applications/CLI/Tests.cmake
+++ b/Applications/CLI/Tests.cmake
@@ -12,7 +12,7 @@ if(NOT OGS_USE_MPI)
EXECUTABLE_ARGS cube_${mesh_size}.prj
WRAPPER time
TESTER vtkdiff
- DIFF_DATA cube_${mesh_size}_pcs_0_ts_1.vtu Linear_1_to_minus1 pressure
+ DIFF_DATA cube_${mesh_size}_pcs_0_ts_1_t_1.000000.vtu Linear_1_to_minus1 pressure
DATA cube_${mesh_size}.prj cube_1x1x1_hex_${mesh_size}.vtu cube_1x1x1.gml
)
@@ -23,7 +23,7 @@ if(NOT OGS_USE_MPI)
EXECUTABLE_ARGS cube_${mesh_size}_neumann.prj
WRAPPER time
TESTER vtkdiff
- DIFF_DATA cube_${mesh_size}_neumann_pcs_0_ts_1.vtu D1_left_front_N1_right pressure
+ DIFF_DATA cube_${mesh_size}_neumann_pcs_0_ts_1_t_1.000000.vtu D1_left_front_N1_right pressure
DATA cube_${mesh_size}_neumann.prj cube_1x1x1_hex_${mesh_size}.vtu cube_1x1x1.gml
)
endforeach()
@@ -36,7 +36,7 @@ if(NOT OGS_USE_MPI)
EXECUTABLE_ARGS cube_${mesh_size}.prj
WRAPPER time
TESTER vtkdiff
- DIFF_DATA cube_${mesh_size}_pcs_0_ts_1.vtu Linear_1_to_minus1 pressure
+ DIFF_DATA cube_${mesh_size}_pcs_0_ts_1_t_1.000000.vtu Linear_1_to_minus1 pressure
DATA cube_${mesh_size}.prj cube_1x1x1_hex_${mesh_size}.vtu cube_1x1x1.gml
)
@@ -47,7 +47,7 @@ if(NOT OGS_USE_MPI)
EXECUTABLE_ARGS cube_${mesh_size}_neumann.prj
WRAPPER time
TESTER vtkdiff
- DIFF_DATA cube_${mesh_size}_neumann_pcs_0_ts_1.vtu D1_left_front_N1_right pressure
+ DIFF_DATA cube_${mesh_size}_neumann_pcs_0_ts_1_t_1.000000.vtu D1_left_front_N1_right pressure
DATA cube_${mesh_size}_neumann.prj cube_1x1x1_hex_${mesh_size}.vtu cube_1x1x1.gml
)
endforeach()
@@ -61,7 +61,7 @@ if(NOT OGS_USE_MPI)
EXECUTABLE_ARGS square_${mesh_size}.prj
WRAPPER time
TESTER vtkdiff
- DIFF_DATA square_${mesh_size}_pcs_0_ts_1.vtu Linear_1_to_minus1 pressure
+ DIFF_DATA square_${mesh_size}_pcs_0_ts_1_t_1.000000.vtu Linear_1_to_minus1 pressure
DATA square_${mesh_size}.prj square_1x1_quad_${mesh_size}.vtu square_1x1.gml
)
@@ -72,7 +72,7 @@ if(NOT OGS_USE_MPI)
EXECUTABLE_ARGS square_${mesh_size}_neumann.prj
WRAPPER time
TESTER vtkdiff
- DIFF_DATA square_${mesh_size}_neumann_pcs_0_ts_1.vtu D1_left_bottom_N1_right pressure
+ DIFF_DATA square_${mesh_size}_neumann_pcs_0_ts_1_t_1.000000.vtu D1_left_bottom_N1_right pressure
DATA square_${mesh_size}_neumann.prj square_1x1_quad_${mesh_size}.vtu square_1x1.gml
)
endforeach()
@@ -85,7 +85,7 @@ if(NOT OGS_USE_MPI)
EXECUTABLE_ARGS square_${mesh_size}.prj
WRAPPER time
TESTER vtkdiff
- DIFF_DATA square_${mesh_size}_pcs_0_ts_1.vtu Linear_1_to_minus1 pressure
+ DIFF_DATA square_${mesh_size}_pcs_0_ts_1_t_1.000000.vtu Linear_1_to_minus1 pressure
DATA square_${mesh_size}.prj square_1x1_quad_${mesh_size}.vtu square_1x1.gml
)
@@ -96,7 +96,7 @@ if(NOT OGS_USE_MPI)
EXECUTABLE_ARGS square_${mesh_size}_neumann.prj
WRAPPER time
TESTER vtkdiff
- DIFF_DATA square_${mesh_size}_neumann_pcs_0_ts_1.vtu D1_left_bottom_N1_right pressure
+ DIFF_DATA square_${mesh_size}_neumann_pcs_0_ts_1_t_1.000000.vtu D1_left_bottom_N1_right pressure
DATA square_${mesh_size}_neumann.prj square_1x1_quad_${mesh_size}.vtu square_1x1.gml
)
endforeach()
@@ -110,7 +110,7 @@ if(NOT OGS_USE_MPI)
EXECUTABLE_ARGS line_${mesh_size}.prj
WRAPPER time
TESTER vtkdiff
- DIFF_DATA line_${mesh_size}_pcs_0_ts_1.vtu Linear_1_to_minus1 pressure
+ DIFF_DATA line_${mesh_size}_pcs_0_ts_1_t_1.000000.vtu Linear_1_to_minus1 pressure
DATA line_${mesh_size}.prj line_1_line_${mesh_size}.vtu line_1.gml
)
@@ -121,7 +121,7 @@ if(NOT OGS_USE_MPI)
EXECUTABLE_ARGS line_${mesh_size}_neumann.prj
WRAPPER time
TESTER vtkdiff
- DIFF_DATA line_${mesh_size}_neumann_pcs_0_ts_1.vtu D1_left_N1_right pressure
+ DIFF_DATA line_${mesh_size}_neumann_pcs_0_ts_1_t_1.000000.vtu D1_left_N1_right pressure
DATA line_${mesh_size}_neumann.prj line_1_line_${mesh_size}.vtu line_1.gml
)
endforeach()
@@ -135,9 +135,9 @@ else()
WRAPPER_ARGS -np 3
TESTER diff
DIFF_DATA
- quad_20x10_GroundWaterFlow_result_pcs_0_ts_1_0.vtu
- quad_20x10_GroundWaterFlow_result_pcs_0_ts_1_1.vtu
- quad_20x10_GroundWaterFlow_result_pcs_0_ts_1_2.vtu
+ quad_20x10_GroundWaterFlow_result_pcs_0_ts_1_t_1.000000_0.vtu
+ quad_20x10_GroundWaterFlow_result_pcs_0_ts_1_t_1.000000_1.vtu
+ quad_20x10_GroundWaterFlow_result_pcs_0_ts_1_t_1.000000_2.vtu
)
AddTest(
@@ -148,9 +148,9 @@ else()
WRAPPER_ARGS -np 3
TESTER diff
DIFF_DATA
- cube_1e3_pcs_0_ts_1_0.vtu
- cube_1e3_pcs_0_ts_1_1.vtu
- cube_1e3_pcs_0_ts_1_2.vtu
+ cube_1e3_pcs_0_ts_1_t_1.000000_0.vtu
+ cube_1e3_pcs_0_ts_1_t_1.000000_1.vtu
+ cube_1e3_pcs_0_ts_1_t_1.000000_2.vtu
)
AddTest(
@@ -161,8 +161,8 @@ else()
WRAPPER_ARGS -np 3
TESTER diff
DIFF_DATA
- cube_1e3_neumann_pcs_0_ts_1_0.vtu
- cube_1e3_neumann_pcs_0_ts_1_1.vtu
- cube_1e3_neumann_pcs_0_ts_1_2.vtu
+ cube_1e3_neumann_pcs_0_ts_1_t_1.000000_0.vtu
+ cube_1e3_neumann_pcs_0_ts_1_t_1.000000_1.vtu
+ cube_1e3_neumann_pcs_0_ts_1_t_1.000000_2.vtu
)
endif()
diff --git a/Applications/CLI/ogs.cpp b/Applications/CLI/ogs.cpp
index d46243b5c68694d85eb2c18e73f1201183e3099f..3702a7ce167240e4e4c54865e5a551bedb039dc9 100644
--- a/Applications/CLI/ogs.cpp
+++ b/Applications/CLI/ogs.cpp
@@ -81,7 +81,8 @@ int main(int argc, char *argv[])
auto project_config = BaseLib::makeConfigTree(
project_arg.getValue(), !nonfatal_arg.getValue(), "OpenGeoSysProject");
- ProjectData project(*project_config, BaseLib::extractPath(project_arg.getValue()));
+ ProjectData project(*project_config, BaseLib::extractPath(project_arg.getValue()),
+ outdir_arg.getValue());
project_config.checkAndInvalidate();
@@ -99,7 +100,7 @@ int main(int argc, char *argv[])
INFO("Solve processes.");
auto& time_loop = project.getTimeLoop();
- time_loop.loop(project, outdir_arg.getValue());
+ time_loop.loop(project);
return 0;
}
diff --git a/AssemblerLib/LocalDataInitializer.h b/AssemblerLib/LocalDataInitializer.h
index 46f69d426960b4194af1976edd1956cff138b14a..7dd42997176d425bbd00968b67c5790cf20c023e 100644
--- a/AssemblerLib/LocalDataInitializer.h
+++ b/AssemblerLib/LocalDataInitializer.h
@@ -111,6 +111,20 @@ static_assert(false, "The macro OGS_MAX_ELEMENT_ORDER is undefined.");
#endif
+// Generates a lambda that creates a new LocalAssembler of type LAData<SHAPE_FCT>
+#define OGS_NEW_LOCAL_ASSEMBLER(SHAPE_FCT) \
+ [](MeshLib::Element const& e, \
+ std::size_t const local_matrix_size, \
+ unsigned const integration_order, \
+ ConstructorArgs&&... args) \
+ { \
+ return new LAData<SHAPE_FCT>{ \
+ e, local_matrix_size, integration_order, \
+ std::forward<ConstructorArgs>(args)... \
+ }; \
+ }
+
+
namespace AssemblerLib
{
@@ -124,7 +138,8 @@ template <
template <typename, typename, typename, typename, unsigned> class LocalAssemblerData_,
typename GlobalMatrix_,
typename GlobalVector_,
- unsigned GlobalDim>
+ unsigned GlobalDim,
+ typename... ConstructorArgs>
class LocalDataInitializer
{
template <typename ShapeFunction_>
@@ -137,7 +152,6 @@ class LocalDataInitializer
IntegrationMethod<ShapeFunction_>,
GlobalMatrix_, GlobalVector_, GlobalDim>;
-
public:
LocalDataInitializer()
{
@@ -146,19 +160,19 @@ public:
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 0 && OGS_MAX_ELEMENT_ORDER >= 1
_builder[std::type_index(typeid(MeshLib::Point))] =
- [](){ return new LAData<NumLib::ShapePoint1>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapePoint1);
#endif
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 1 && OGS_MAX_ELEMENT_ORDER >= 1
_builder[std::type_index(typeid(MeshLib::Line))] =
- [](){ return new LAData<NumLib::ShapeLine2>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapeLine2);
#endif
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 1 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Line3))] =
- [](){ return new LAData<NumLib::ShapeLine3>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapeLine3);
#endif
@@ -167,27 +181,27 @@ public:
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 1
_builder[std::type_index(typeid(MeshLib::Quad))] =
- [](){ return new LAData<NumLib::ShapeQuad4>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapeQuad4);
#endif
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1
_builder[std::type_index(typeid(MeshLib::Hex))] =
- [](){ return new LAData<NumLib::ShapeHex8>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapeHex8);
#endif
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Quad8))] =
- [](){ return new LAData<NumLib::ShapeQuad8>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapeQuad8);
_builder[std::type_index(typeid(MeshLib::Quad9))] =
- [](){ return new LAData<NumLib::ShapeQuad9>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapeQuad9);
#endif
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Hex20))] =
- [](){ return new LAData<NumLib::ShapeHex20>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapeHex20);
#endif
@@ -196,25 +210,25 @@ public:
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 1
_builder[std::type_index(typeid(MeshLib::Tri))] =
- [](){ return new LAData<NumLib::ShapeTri3>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapeTri3);
#endif
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1
_builder[std::type_index(typeid(MeshLib::Tet))] =
- [](){ return new LAData<NumLib::ShapeTet4>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapeTet4);
#endif
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Tri6))] =
- [](){ return new LAData<NumLib::ShapeTri6>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapeTri6);
#endif
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Tet10))] =
- [](){ return new LAData<NumLib::ShapeTet10>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapeTet10);
#endif
@@ -223,13 +237,13 @@ public:
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1
_builder[std::type_index(typeid(MeshLib::Prism))] =
- [](){ return new LAData<NumLib::ShapePrism6>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapePrism6);
#endif
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Prism15))] =
- [](){ return new LAData<NumLib::ShapePrism15>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapePrism15);
#endif
// /// Pyramids //////////////////////////////////////////////////
@@ -237,28 +251,31 @@ public:
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1
_builder[std::type_index(typeid(MeshLib::Pyramid))] =
- [](){ return new LAData<NumLib::ShapePyra5>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapePyra5);
#endif
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Pyramid13))] =
- [](){ return new LAData<NumLib::ShapePyra13>; };
+ OGS_NEW_LOCAL_ASSEMBLER(NumLib::ShapePyra13);
#endif
}
- /// Sets the provided data_ptr to the newly created local assembler data and
- /// calls init() forwarding all remaining arguments.
- template <typename ...Args_>
- void operator()(const MeshLib::Element& e,
- LocalAssemblerDataInterface_<GlobalMatrix_, GlobalVector_>*& data_ptr, Args_&&... args)
+ /// Sets the provided data_ptr to the newly created local assembler data.
+ void operator()(
+ const MeshLib::Element& e,
+ LocalAssemblerDataInterface_<GlobalMatrix_, GlobalVector_>*& data_ptr,
+ std::size_t const local_matrix_size,
+ unsigned const integration_order,
+ ConstructorArgs&&... args)
{
auto const type_idx = std::type_index(typeid(e));
auto it = _builder.find(type_idx);
if (it != _builder.end()) {
- data_ptr = it->second();
- data_ptr->init(e, std::forward<Args_>(args)...);
+ data_ptr = it->second(
+ e, local_matrix_size, integration_order,
+ std::forward<ConstructorArgs>(args)...);
} else {
ERR("You are trying to build a local assembler for an unknown mesh element type (%s)."
" Maybe you have disabled this mesh element type in your build configuration.",
@@ -271,8 +288,12 @@ private:
/// Mapping of element types to local assembler constructors.
std::unordered_map<
std::type_index,
- std::function<LocalAssemblerDataInterface_<GlobalMatrix_, GlobalVector_>*()>
- > _builder;
+ std::function<LocalAssemblerDataInterface_<GlobalMatrix_, GlobalVector_>*(
+ MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ unsigned const integration_order,
+ ConstructorArgs&&...)>
+ > _builder;
};
} // namespace AssemblerLib
@@ -286,5 +307,6 @@ private:
#undef ENABLED_ELEMENT_TYPE_TRI
#undef ENABLED_ELEMENT_TYPE_QUAD
#undef OGS_ENABLED_ELEMENTS
+#undef OGS_NEW_LOCAL_ASSEMBLER
#endif // ASSEMBLER_LIB_LOCALDATAINITIALIZER_H_
diff --git a/FileIO/VtkIO/PVDFile.cpp b/FileIO/VtkIO/PVDFile.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..14a6430aa237d13456d8fede8afbd6efec76600c
--- /dev/null
+++ b/FileIO/VtkIO/PVDFile.cpp
@@ -0,0 +1,43 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "PVDFile.h"
+
+#include <fstream>
+#include <iomanip>
+#include <limits>
+#include <logog/include/logog.hpp>
+
+namespace FileIO
+{
+
+void PVDFile::addVTUFile(const std::string &vtu_fname, double timestep)
+{
+ _datasets.push_back(std::make_pair(timestep, vtu_fname));
+
+ std::ofstream fh(_pvd_filename.c_str());
+ if (!fh) {
+ ERR("could not open file `%s'", _pvd_filename.c_str());
+ std::abort();
+ }
+
+ fh << std::setprecision(std::numeric_limits<double>::digits10);
+
+ fh << "<?xml version=\"1.0\"?>\n"
+ "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\""
+ " compressor=\"vtkZLibDataCompressor\">\n"
+ " <Collection>\n";
+
+ for (auto const& pair : _datasets)
+ fh << " <DataSet timestep=\"" << pair.first << "\" group=\"\" part=\"0\" file=\"" << pair.second << "\"/>\n";
+
+ fh << " </Collection>\n</VTKFile>\n";
+}
+
+}
diff --git a/FileIO/VtkIO/PVDFile.h b/FileIO/VtkIO/PVDFile.h
new file mode 100644
index 0000000000000000000000000000000000000000..4cf48618d62daa5a77b8102c3c8dd3e9cb4d91ec
--- /dev/null
+++ b/FileIO/VtkIO/PVDFile.h
@@ -0,0 +1,38 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef FILEIO_VTK_PVDFILE_H
+#define FILEIO_VTK_PVDFILE_H
+
+#include <string>
+#include <vector>
+
+namespace FileIO
+{
+
+/*! Writes a basic PVD file for use with Paraview.
+ *
+ */
+class PVDFile
+{
+public:
+ //! Set a PVD file path
+ explicit PVDFile(std::string const& pvd_fname) : _pvd_filename(pvd_fname) {}
+
+ //! Add a VTU file to this PVD file.
+ void addVTUFile(std::string const& vtu_fname, double timestep);
+
+private:
+ std::string const _pvd_filename;
+ std::vector<std::pair<double, std::string>> _datasets; // a vector of (time, VTU file name)
+};
+
+} // namespace FileIO
+
+#endif // FILEIO_VTK_PVDFILE_H
diff --git a/ProcessLib/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlowFEM.h
index 6fabf1b3ea9929c6d57ceafc424158f615b9e965..812c0496be55632d33ce9b99b90a270cfc1e3450 100644
--- a/ProcessLib/GroundwaterFlowFEM.h
+++ b/ProcessLib/GroundwaterFlowFEM.h
@@ -32,11 +32,6 @@ class LocalAssemblerDataInterface
public:
virtual ~LocalAssemblerDataInterface() = default;
- virtual void init(MeshLib::Element const& e,
- std::size_t const local_matrix_size,
- Parameter<double, MeshLib::Element const&> const& hydraulic_conductivity,
- unsigned const integration_order) = 0;
-
virtual void assemble(double const t, std::vector<double> const& local_x) = 0;
virtual void addToGlobal(AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const&,
@@ -59,11 +54,11 @@ public:
/// The hydraulic_conductivity factor is directly integrated into the local
/// element matrix.
- void init(MeshLib::Element const& e,
- std::size_t const local_matrix_size,
- Parameter<double, MeshLib::Element const&> const&
- hydraulic_conductivity,
- unsigned const integration_order) override
+ LocalAssemblerData(MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ unsigned const integration_order,
+ Parameter<double, MeshLib::Element const&> const&
+ hydraulic_conductivity)
{
_integration_order = integration_order;
diff --git a/ProcessLib/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlowProcess.h
index e7641fafcec211305eb7d59c81a95369af6fb527..528711d1f879e6be162324a3f97ef8fc4b7f2e87 100644
--- a/ProcessLib/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlowProcess.h
@@ -80,7 +80,8 @@ public:
GroundwaterFlow::LocalAssemblerData,
typename GlobalSetup::MatrixType,
typename GlobalSetup::VectorType,
- GlobalDim>;
+ GlobalDim,
+ Parameter<double, MeshLib::Element const&> const&>;
LocalDataInitializer initializer;
@@ -97,8 +98,8 @@ public:
local_asm_builder,
this->_mesh.getElements(),
_local_assemblers,
- _hydraulic_conductivity,
- this->_integration_order);
+ this->_integration_order,
+ _hydraulic_conductivity);
}
// TODO remove, but put "gw_" somewhere
diff --git a/ProcessLib/NeumannBc.h b/ProcessLib/NeumannBc.h
index d9faba93024aa1aa0138891d2d407bd658c29426..8cc1349e0edd4551014c1de20f33a4cdf4bc969f 100644
--- a/ProcessLib/NeumannBc.h
+++ b/ProcessLib/NeumannBc.h
@@ -131,7 +131,8 @@ private:
LocalNeumannBcAsmDataInterface,
LocalNeumannBcAsmData,
GlobalMatrix, GlobalVector,
- GlobalDim>;
+ GlobalDim,
+ std::function<double (MeshLib::Element const&)> const&>;
LocalDataInitializer initializer;
@@ -155,8 +156,8 @@ private:
local_asm_builder,
_elements,
_local_assemblers,
- elementValueLookup,
- _integration_order);
+ _integration_order,
+ elementValueLookup);
DBUG("Create global assembler.");
_global_assembler.reset(
diff --git a/ProcessLib/NeumannBcAssembler.h b/ProcessLib/NeumannBcAssembler.h
index 9e610e15ce3b31378de7aabef01e38fefa7b857f..fe124ef9996aa3659552fd23af1c3e0eddf672c1 100644
--- a/ProcessLib/NeumannBcAssembler.h
+++ b/ProcessLib/NeumannBcAssembler.h
@@ -25,11 +25,6 @@ class LocalNeumannBcAsmDataInterface
public:
virtual ~LocalNeumannBcAsmDataInterface() = default;
- virtual void init(MeshLib::Element const& e,
- std::size_t const local_matrix_size,
- std::function<double (MeshLib::Element const&)> const& value_lookup,
- unsigned const integration_order) = 0;
-
virtual void assemble(const double t) = 0;
virtual void addToGlobal(AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const&,
@@ -53,11 +48,11 @@ public:
/// The neumann_bc_value factor is directly integrated into the local
/// element matrix.
- void
- init(MeshLib::Element const& e,
- std::size_t const local_matrix_size,
- std::function<double (MeshLib::Element const&)> const& value_lookup,
- unsigned const integration_order) override
+ LocalNeumannBcAsmData(
+ MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ unsigned const integration_order,
+ std::function<double (MeshLib::Element const&)> const& value_lookup)
{
using FemType = NumLib::TemplateIsoparametric<
ShapeFunction, ShapeMatricesType>;
diff --git a/ProcessLib/Output.cpp b/ProcessLib/Output.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0646f2bd9a5fdd07b19d4473ed0fcd28f33a3744
--- /dev/null
+++ b/ProcessLib/Output.cpp
@@ -0,0 +1,147 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include<fstream>
+#include "Output.h"
+
+#include<cassert>
+#include<vector>
+
+#include <logog/include/logog.hpp>
+
+#include "BaseLib/FileTools.h"
+
+namespace
+{
+//! Determines if there should be output at the given \c timestep.
+template<typename CountsSteps>
+bool shallDoOutput(unsigned timestep, CountsSteps const& repeats_each_steps)
+{
+ unsigned each_steps = 1;
+
+ for (auto const& pair : repeats_each_steps)
+ {
+ each_steps = pair.each_steps;
+
+ if (timestep > pair.repeat * each_steps)
+ {
+ timestep -= pair.repeat * each_steps;
+ } else {
+ break;
+ }
+ }
+
+ return timestep % each_steps == 0;
+}
+}
+
+namespace ProcessLib
+{
+
+template<typename GlobalSetup>
+std::unique_ptr<Output<GlobalSetup>> Output<GlobalSetup>::
+newInstance(const BaseLib::ConfigTree &config, std::string const& output_directory)
+{
+ std::unique_ptr<Output> out{ new Output{
+ BaseLib::joinPaths(output_directory,
+ config.getConfParam<std::string>("prefix"))}};
+
+ if (auto const timesteps = config.getConfSubtreeOptional("timesteps"))
+ {
+ for (auto pair : timesteps->getConfSubtreeList("pair"))
+ {
+ auto repeat = pair.getConfParam<unsigned>("repeat");
+ auto each_steps = pair.getConfParam<unsigned>("each_steps");
+
+ assert(repeat != 0 && each_steps != 0);
+ out->_repeats_each_steps.emplace_back(repeat, each_steps);
+ }
+
+ if (out->_repeats_each_steps.empty()) {
+ ERR("You have not given any pair (<repeat/>, <each_steps/>) that defines"
+ " at which timesteps output shall be written. Aborting.");
+ std::abort();
+ }
+ }
+ else
+ {
+ out->_repeats_each_steps.emplace_back(1, 1);
+ }
+
+ return out;
+}
+
+template<typename GlobalSetup>
+void Output<GlobalSetup>::
+init(typename Output<GlobalSetup>::ProcessIter first,
+ typename Output<GlobalSetup>::ProcessIter last)
+{
+ for (unsigned pcs_idx = 0; first != last; ++first, ++pcs_idx)
+ {
+ auto const filename = _output_file_prefix
+ + "_pcs_" + std::to_string(pcs_idx)
+ + ".pvd";
+ _single_process_data.emplace(std::piecewise_construct,
+ std::forward_as_tuple(&**first),
+ std::forward_as_tuple(pcs_idx, filename));
+ }
+}
+
+template<typename GlobalSetup>
+void Output<GlobalSetup>::
+doOutputAlways(Process<GlobalSetup> const& process,
+ unsigned timestep,
+ const double t,
+ typename GlobalSetup::VectorType const& x)
+{
+ auto spd_it = _single_process_data.find(&process);
+ if (spd_it == _single_process_data.end()) {
+ ERR("The given process is not contained in the output configuration."
+ " Aborting.");
+ std::abort();
+ }
+ auto& spd = spd_it->second;
+
+ std::string const output_file_name =
+ _output_file_prefix + "_pcs_" + std::to_string(spd.process_index)
+ + "_ts_" + std::to_string(timestep)
+ + "_t_" + std::to_string(t)
+ + ".vtu";
+ DBUG("output to %s", output_file_name.c_str());
+ process.output(output_file_name, timestep, x);
+ spd.pvd_file.addVTUFile(output_file_name, t);
+}
+
+template<typename GlobalSetup>
+void Output<GlobalSetup>::
+doOutput(Process<GlobalSetup> const& process,
+ unsigned timestep,
+ const double t,
+ typename GlobalSetup::VectorType const& x)
+{
+ if (shallDoOutput(timestep, _repeats_each_steps))
+ doOutputAlways(process, timestep, t, x);
+}
+
+template<typename GlobalSetup>
+void Output<GlobalSetup>::
+doOutputLastTimestep(Process<GlobalSetup> const& process,
+ unsigned timestep,
+ const double t,
+ typename GlobalSetup::VectorType const& x)
+{
+ if (!shallDoOutput(timestep, _repeats_each_steps))
+ doOutputAlways(process, timestep, t, x);
+}
+
+
+// explicit instantiation
+template class Output<GlobalSetupType>;
+
+}
diff --git a/ProcessLib/Output.h b/ProcessLib/Output.h
new file mode 100644
index 0000000000000000000000000000000000000000..b83826acac758cca7a43419a5ff933f6f0eb6a43
--- /dev/null
+++ b/ProcessLib/Output.h
@@ -0,0 +1,98 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_OUTPUT_H
+#define PROCESSLIB_OUTPUT_H
+
+#include "BaseLib/ConfigTree.h"
+#include "FileIO/VtkIO/PVDFile.h"
+#include "Process.h"
+
+namespace ProcessLib
+{
+
+/*! Manages writing the solution of processes to disk.
+ *
+ * This class decides at which timesteps output is written
+ * and initiates the writing process.
+ */
+template<typename GlobalSetup>
+class Output
+{
+public:
+ static std::unique_ptr<Output>
+ newInstance(const BaseLib::ConfigTree& config,
+ const std::string& output_directory);
+
+ using ProcessIter = std::vector<std::unique_ptr<ProcessLib::Process<GlobalSetupType>>>
+ ::const_iterator;
+
+ //! Opens a PVD file for each process.
+ void init(ProcessIter first, ProcessIter last);
+
+ //! Writes output for the given \c process if it should be written in the
+ //! given \c timestep.
+ void doOutput(
+ Process<GlobalSetup> const& process, unsigned timestep,
+ const double t,
+ typename GlobalSetup::VectorType const& x);
+
+ //! Writes output for the given \c process if it has not been written yet.
+ //! This method is intended for doing output after the last timestep in order
+ //! to make sure that its results are written.
+ void doOutputLastTimestep(
+ Process<GlobalSetup> const& process, unsigned timestep,
+ const double t,
+ typename GlobalSetup::VectorType const& x);
+
+ //! Writes output for the given \c process.
+ //! This method will always write.
+ //! It is intended to write output in error handling routines.
+ void doOutputAlways(
+ Process<GlobalSetup> const& process, unsigned timestep,
+ const double t,
+ typename GlobalSetup::VectorType const& x);
+
+ struct PairRepeatEachSteps
+ {
+ explicit PairRepeatEachSteps(unsigned c, unsigned e)
+ : repeat(c), each_steps(e)
+ {}
+
+ const unsigned repeat; //!< Apply \c each_steps \c repeat times.
+ const unsigned each_steps; //!< Do output every \c each_steps timestep.
+ };
+private:
+ struct SingleProcessData
+ {
+ SingleProcessData(unsigned process_index_,
+ std::string const& filename)
+ : process_index(process_index_)
+ , pvd_file(filename)
+ {}
+
+ const unsigned process_index;
+ FileIO::PVDFile pvd_file;
+ };
+
+ Output(std::string const& prefix)
+ : _output_file_prefix(prefix)
+ {}
+
+ std::string _output_file_prefix;
+
+ //! Describes after which timesteps to write output.
+ std::vector<PairRepeatEachSteps> _repeats_each_steps;
+
+ std::map<Process<GlobalSetup> const*, SingleProcessData> _single_process_data;
+};
+
+}
+
+#endif // PROCESSLIB_OUTPUT_H
diff --git a/Tests/Data b/Tests/Data
index 879c870ddea452ecf3d8e0bd5a769d15b73f8be1..2464c6e74e85924556091735120db4cedf6abdf6 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit 879c870ddea452ecf3d8e0bd5a769d15b73f8be1
+Subproject commit 2464c6e74e85924556091735120db4cedf6abdf6
diff --git a/scripts/cmake/test/AddTest.cmake b/scripts/cmake/test/AddTest.cmake
index db228b0d42b66c2c194bc50379720eb840711d56..66a9aa3175a085c0a9479b7a3dbab05e8ccac87e 100644
--- a/scripts/cmake/test/AddTest.cmake
+++ b/scripts/cmake/test/AddTest.cmake
@@ -19,11 +19,11 @@
#
# diff-tester
# - DIFF_DATA <list of files to diff>
-# # the given file is compared to [filename]_expected.[extension]
+# # the given file is compared to a file with the same name from Tests/Data
#
# numdiff-tester
# - DIFF_DATA <list of files to numdiff>
-# # the given file is compared to [filename]_expected.[extension]
+# # the given file is compared to a file with the same name from Tests/Data
#
# vtkdiff-tester
# - DIFF_DATA <vtk file> <data array a name> <data array b name>
@@ -115,9 +115,7 @@ function (AddTest)
if(AddTest_TESTER STREQUAL "diff" OR AddTest_TESTER STREQUAL "numdiff")
foreach(FILE ${AddTest_DIFF_DATA})
- get_filename_component(FILE_NAME ${FILE} NAME_WE)
- get_filename_component(FILE_EXT ${FILE} EXT)
- set(FILE_EXPECTED ${FILE_NAME}_expected${FILE_EXT})
+ get_filename_component(FILE_EXPECTED ${FILE} NAME)
set(TESTER_COMMAND ${TESTER_COMMAND} "${SELECTED_DIFF_TOOL_PATH} \
${TESTER_ARGS} ${AddTest_SOURCE_PATH}/${FILE_EXPECTED} \
${AddTest_BINARY_PATH}/${FILE}")