Skip to content
Snippets Groups Projects
Commit 7bb833c4 authored by Dmitri Naumov's avatar Dmitri Naumov
Browse files

[PL] Tabs to whitespaces.

parent 48551248
No related branches found
No related tags found
No related merge requests found
...@@ -18,8 +18,8 @@ target_link_libraries(ProcessLib ...@@ -18,8 +18,8 @@ target_link_libraries(ProcessLib
ADD_VTK_DEPENDENCY(ProcessLib) ADD_VTK_DEPENDENCY(ProcessLib)
if(TARGET Eigen) if(TARGET Eigen)
add_dependencies(ProcessLib Eigen) add_dependencies(ProcessLib Eigen)
endif() endif()
if(TARGET Boost) if(TARGET Boost)
add_dependencies(ProcessLib Boost) add_dependencies(ProcessLib Boost)
endif() endif()
...@@ -22,8 +22,8 @@ namespace ProcessLib ...@@ -22,8 +22,8 @@ namespace ProcessLib
template <typename IndexType> template <typename IndexType>
struct DirichletBc final struct DirichletBc final
{ {
std::vector<IndexType> global_ids; std::vector<IndexType> global_ids;
std::vector<double> values; std::vector<double> values;
}; };
} // namespace ProcessLib } // namespace ProcessLib
......
...@@ -23,13 +23,13 @@ namespace ProcessLib ...@@ -23,13 +23,13 @@ namespace ProcessLib
std::unique_ptr<InitialCondition> createUniformInitialCondition( std::unique_ptr<InitialCondition> createUniformInitialCondition(
BaseLib::ConfigTree const& config, int const /*n_components*/) BaseLib::ConfigTree const& config, int const /*n_components*/)
{ {
config.checkConfParam("type", "Uniform"); config.checkConfParam("type", "Uniform");
auto value = config.getConfParam<double>("value"); auto value = config.getConfParam<double>("value");
DBUG("Using value %g", value); DBUG("Using value %g", value);
return std::unique_ptr<InitialCondition>( return std::unique_ptr<InitialCondition>(
new UniformInitialCondition(value)); new UniformInitialCondition(value));
} }
std::unique_ptr<InitialCondition> createMeshPropertyInitialCondition( std::unique_ptr<InitialCondition> createMeshPropertyInitialCondition(
...@@ -37,34 +37,34 @@ std::unique_ptr<InitialCondition> createMeshPropertyInitialCondition( ...@@ -37,34 +37,34 @@ std::unique_ptr<InitialCondition> createMeshPropertyInitialCondition(
MeshLib::Mesh const& mesh, MeshLib::Mesh const& mesh,
int const n_components) int const n_components)
{ {
auto field_name = config.getConfParam<std::string>("field_name"); auto field_name = config.getConfParam<std::string>("field_name");
DBUG("Using field_name %s", field_name.c_str()); DBUG("Using field_name %s", field_name.c_str());
if (!mesh.getProperties().hasPropertyVector(field_name)) if (!mesh.getProperties().hasPropertyVector(field_name))
{ {
ERR("The required property %s does not exists in the mesh.", ERR("The required property %s does not exists in the mesh.",
field_name.c_str()); field_name.c_str());
std::abort(); std::abort();
} }
auto const& property = auto const& property =
mesh.getProperties().template getPropertyVector<double>(field_name); mesh.getProperties().template getPropertyVector<double>(field_name);
if (!property) if (!property)
{ {
ERR("The required property %s is not of the requested type.", ERR("The required property %s is not of the requested type.",
field_name.c_str()); field_name.c_str());
std::abort(); std::abort();
} }
if (property->getNumberOfComponents() != if (property->getNumberOfComponents() !=
static_cast<std::size_t>(n_components)) static_cast<std::size_t>(n_components))
{ {
ERR("The required property %s has different number of components %d, " ERR("The required property %s has different number of components %d, "
"expected %d.", "expected %d.",
field_name.c_str(), property->getNumberOfComponents(), n_components); field_name.c_str(), property->getNumberOfComponents(), n_components);
std::abort(); std::abort();
} }
return std::unique_ptr<InitialCondition>( return std::unique_ptr<InitialCondition>(
new MeshPropertyInitialCondition(*property)); new MeshPropertyInitialCondition(*property));
} }
} // namespace ProcessLib } // namespace ProcessLib
...@@ -34,27 +34,27 @@ namespace ProcessLib ...@@ -34,27 +34,27 @@ namespace ProcessLib
class InitialCondition class InitialCondition
{ {
public: public:
virtual ~InitialCondition() = default; virtual ~InitialCondition() = default;
virtual double getValue(MeshLib::Node const&, int const component_id) const = 0; virtual double getValue(MeshLib::Node const&, int const component_id) const = 0;
}; };
/// Uniform value initial condition /// Uniform value initial condition
class UniformInitialCondition : public InitialCondition class UniformInitialCondition : public InitialCondition
{ {
public: public:
UniformInitialCondition(double const value) : _value(value) UniformInitialCondition(double const value) : _value(value)
{ {
} }
/// Returns a value for given node and component. /// Returns a value for given node and component.
/// \todo The component_id is to be implemented. /// \todo The component_id is to be implemented.
virtual double getValue(MeshLib::Node const&, virtual double getValue(MeshLib::Node const&,
int const /* component_id */) const override int const /* component_id */) const override
{ {
return _value; return _value;
} }
private: private:
double _value; double _value;
}; };
/// Construct a UniformInitialCondition from configuration. /// Construct a UniformInitialCondition from configuration.
...@@ -68,22 +68,22 @@ std::unique_ptr<InitialCondition> createUniformInitialCondition( ...@@ -68,22 +68,22 @@ std::unique_ptr<InitialCondition> createUniformInitialCondition(
class MeshPropertyInitialCondition : public InitialCondition class MeshPropertyInitialCondition : public InitialCondition
{ {
public: public:
MeshPropertyInitialCondition( MeshPropertyInitialCondition(
MeshLib::PropertyVector<double> const& property) MeshLib::PropertyVector<double> const& property)
: _property(property) : _property(property)
{ {
assert(_property.getMeshItemType() == MeshLib::MeshItemType::Node); assert(_property.getMeshItemType() == MeshLib::MeshItemType::Node);
} }
virtual double getValue(MeshLib::Node const& n, virtual double getValue(MeshLib::Node const& n,
int const component_id) const override int const component_id) const override
{ {
return _property[n.getID() * _property.getNumberOfComponents() + return _property[n.getID() * _property.getNumberOfComponents() +
component_id]; component_id];
} }
private: private:
MeshLib::PropertyVector<double> const& _property; MeshLib::PropertyVector<double> const& _property;
}; };
/// Construct a MeshPropertyInitialCondition from configuration. /// Construct a MeshPropertyInitialCondition from configuration.
......
...@@ -19,36 +19,36 @@ namespace ProcessLib ...@@ -19,36 +19,36 @@ namespace ProcessLib
std::unique_ptr<ParameterBase> createConstParameter( std::unique_ptr<ParameterBase> createConstParameter(
BaseLib::ConfigTree const& config) BaseLib::ConfigTree const& config)
{ {
config.checkConfParam("type", "Constant"); config.checkConfParam("type", "Constant");
auto value = config.getConfParam<double>("value"); auto value = config.getConfParam<double>("value");
DBUG("Using value %g", value); DBUG("Using value %g", value);
return std::unique_ptr<ParameterBase>(new ConstParameter<double>(value)); return std::unique_ptr<ParameterBase>(new ConstParameter<double>(value));
} }
std::unique_ptr<ParameterBase> createMeshPropertyParameter( std::unique_ptr<ParameterBase> createMeshPropertyParameter(
BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh) BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh)
{ {
config.checkConfParam("type", "MeshProperty"); config.checkConfParam("type", "MeshProperty");
auto field_name = config.getConfParam<std::string>("field_name"); auto field_name = config.getConfParam<std::string>("field_name");
DBUG("Using field_name %s", field_name.c_str()); DBUG("Using field_name %s", field_name.c_str());
if (!mesh.getProperties().hasPropertyVector(field_name)) if (!mesh.getProperties().hasPropertyVector(field_name))
{ {
ERR("The required property %s does not exists in the mesh.", ERR("The required property %s does not exists in the mesh.",
field_name.c_str()); field_name.c_str());
std::abort(); std::abort();
} }
auto const& property = auto const& property =
mesh.getProperties().template getPropertyVector<double>(field_name); mesh.getProperties().template getPropertyVector<double>(field_name);
if (!property) if (!property)
{ {
ERR("The required property %s is not of the requested type.", ERR("The required property %s is not of the requested type.",
field_name.c_str()); field_name.c_str());
std::abort(); std::abort();
} }
return std::unique_ptr<ParameterBase>( return std::unique_ptr<ParameterBase>(
new MeshPropertyParameter<double>(*property)); new MeshPropertyParameter<double>(*property));
} }
} // namespace ProcessLib } // namespace ProcessLib
...@@ -31,9 +31,9 @@ namespace ProcessLib ...@@ -31,9 +31,9 @@ namespace ProcessLib
/// Its property name helps addressing the right parameter. /// Its property name helps addressing the right parameter.
struct ParameterBase struct ParameterBase
{ {
virtual ~ParameterBase() = default; virtual ~ParameterBase() = default;
std::string name; std::string name;
}; };
/// A parameter is representing a value or function of any type. /// A parameter is representing a value or function of any type.
...@@ -43,9 +43,9 @@ struct ParameterBase ...@@ -43,9 +43,9 @@ struct ParameterBase
template <typename ReturnType, typename... Args> template <typename ReturnType, typename... Args>
struct Parameter : public ParameterBase struct Parameter : public ParameterBase
{ {
virtual ~Parameter() = default; virtual ~Parameter() = default;
virtual ReturnType operator()(Args&&... args) const = 0; virtual ReturnType operator()(Args&&... args) const = 0;
}; };
/// Single, constant value parameter. /// Single, constant value parameter.
...@@ -53,17 +53,17 @@ template <typename ReturnType> ...@@ -53,17 +53,17 @@ template <typename ReturnType>
struct ConstParameter final struct ConstParameter final
: public Parameter<ReturnType, MeshLib::Element const&> : public Parameter<ReturnType, MeshLib::Element const&>
{ {
ConstParameter(ReturnType value) : _value(value) ConstParameter(ReturnType value) : _value(value)
{ {
} }
ReturnType operator()(MeshLib::Element const&) const override ReturnType operator()(MeshLib::Element const&) const override
{ {
return _value; return _value;
} }
private: private:
ReturnType _value; ReturnType _value;
}; };
std::unique_ptr<ParameterBase> createConstParameter(BaseLib::ConfigTree const& config); std::unique_ptr<ParameterBase> createConstParameter(BaseLib::ConfigTree const& config);
...@@ -73,18 +73,18 @@ template <typename ReturnType> ...@@ -73,18 +73,18 @@ template <typename ReturnType>
struct MeshPropertyParameter final struct MeshPropertyParameter final
: public Parameter<ReturnType, MeshLib::Element const&> : public Parameter<ReturnType, MeshLib::Element const&>
{ {
MeshPropertyParameter(MeshLib::PropertyVector<ReturnType> const& property) MeshPropertyParameter(MeshLib::PropertyVector<ReturnType> const& property)
: _property(property) : _property(property)
{ {
} }
ReturnType operator()(MeshLib::Element const& e) const override ReturnType operator()(MeshLib::Element const& e) const override
{ {
return _property[e.getID()]; return _property[e.getID()];
} }
private: private:
MeshLib::PropertyVector<ReturnType> const& _property; MeshLib::PropertyVector<ReturnType> const& _property;
}; };
std::unique_ptr<ParameterBase> createMeshPropertyParameter(BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh); std::unique_ptr<ParameterBase> createMeshPropertyParameter(BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh);
......
...@@ -15,47 +15,47 @@ ProcessVariable& findProcessVariable( ...@@ -15,47 +15,47 @@ ProcessVariable& findProcessVariable(
std::vector<ProcessVariable> const& variables, std::vector<ProcessVariable> const& variables,
BaseLib::ConfigTree const& pv_config, std::string const& tag) BaseLib::ConfigTree const& pv_config, std::string const& tag)
{ {
// Find process variable name in process config. // Find process variable name in process config.
std::string const name = pv_config.getConfParam<std::string>(tag); std::string const name = pv_config.getConfParam<std::string>(tag);
// Find corresponding variable by name. // Find corresponding variable by name.
auto variable = std::find_if(variables.cbegin(), variables.cend(), auto variable = std::find_if(variables.cbegin(), variables.cend(),
[&name](ProcessVariable const& v) [&name](ProcessVariable const& v)
{ {
return v.getName() == name; return v.getName() == name;
}); });
if (variable == variables.end()) if (variable == variables.end())
{ {
ERR( ERR(
"Could not find process variable '%s' in the provided variables " "Could not find process variable '%s' in the provided variables "
"list for config tag <%s>.", "list for config tag <%s>.",
name.c_str(), tag.c_str()); name.c_str(), tag.c_str());
std::abort(); std::abort();
} }
DBUG("Found process variable \'%s\' for config tag <%s>.", DBUG("Found process variable \'%s\' for config tag <%s>.",
variable->getName().c_str(), tag.c_str()); variable->getName().c_str(), tag.c_str());
// Const cast is needed because of variables argument constness. // Const cast is needed because of variables argument constness.
return const_cast<ProcessVariable&>(*variable); return const_cast<ProcessVariable&>(*variable);
} }
std::vector<std::reference_wrapper<ProcessVariable>> std::vector<std::reference_wrapper<ProcessVariable>>
findProcessVariables( findProcessVariables(
std::vector<ProcessVariable> const& variables, std::vector<ProcessVariable> const& variables,
BaseLib::ConfigTree const& process_config, BaseLib::ConfigTree const& process_config,
std::initializer_list<std::string> tag_names) std::initializer_list<std::string> tag_names)
{ {
std::vector<std::reference_wrapper<ProcessVariable>> vars; std::vector<std::reference_wrapper<ProcessVariable>> vars;
vars.reserve(tag_names.size()); vars.reserve(tag_names.size());
auto const pv_conf = process_config.getConfSubtree("process_variables"); auto const pv_conf = process_config.getConfSubtree("process_variables");
for (auto const& tag : tag_names) { for (auto const& tag : tag_names) {
vars.emplace_back(findProcessVariable(variables, pv_conf, tag)); vars.emplace_back(findProcessVariable(variables, pv_conf, tag));
} }
return vars; return vars;
} }
} // namespace ProcessLib } // namespace ProcessLib
...@@ -34,294 +34,294 @@ namespace ProcessLib ...@@ -34,294 +34,294 @@ namespace ProcessLib
template <typename GlobalSetup> template <typename GlobalSetup>
class Process class Process
: public NumLib::ODESystem<typename GlobalSetup::MatrixType, : public NumLib::ODESystem<typename GlobalSetup::MatrixType,
typename GlobalSetup::VectorType, typename GlobalSetup::VectorType,
// TODO: later on use a simpler ODE system // TODO: later on use a simpler ODE system
NumLib::ODESystemTag::FirstOrderImplicitQuasilinear, NumLib::ODESystemTag::FirstOrderImplicitQuasilinear,
NumLib::NonlinearSolverTag::Newton> NumLib::NonlinearSolverTag::Newton>
{ {
public: public:
using GlobalVector = typename GlobalSetup::VectorType; using GlobalVector = typename GlobalSetup::VectorType;
using GlobalMatrix = typename GlobalSetup::MatrixType; using GlobalMatrix = typename GlobalSetup::MatrixType;
using Index = typename GlobalMatrix::IndexType; using Index = typename GlobalMatrix::IndexType;
using NonlinearSolver = NumLib::NonlinearSolverBase<GlobalMatrix, GlobalVector>; using NonlinearSolver = NumLib::NonlinearSolverBase<GlobalMatrix, GlobalVector>;
using TimeDiscretization = NumLib::TimeDiscretization<GlobalVector>; using TimeDiscretization = NumLib::TimeDiscretization<GlobalVector>;
Process( Process(
MeshLib::Mesh& mesh, MeshLib::Mesh& mesh,
NonlinearSolver& nonlinear_solver, NonlinearSolver& nonlinear_solver,
std::unique_ptr<TimeDiscretization>&& time_discretization, std::unique_ptr<TimeDiscretization>&& time_discretization,
std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables, std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
SecondaryVariableCollection<GlobalVector>&& secondary_variables, SecondaryVariableCollection<GlobalVector>&& secondary_variables,
ProcessOutput<GlobalVector>&& process_output ProcessOutput<GlobalVector>&& process_output
) )
: _mesh(mesh) : _mesh(mesh)
, _nonlinear_solver(nonlinear_solver) , _nonlinear_solver(nonlinear_solver)
, _time_discretization(std::move(time_discretization)) , _time_discretization(std::move(time_discretization))
, _process_variables(std::move(process_variables)) , _process_variables(std::move(process_variables))
, _process_output(std::move(process_output)) , _process_output(std::move(process_output))
, _secondary_variables(std::move(secondary_variables)) , _secondary_variables(std::move(secondary_variables))
{} {}
/// Preprocessing before starting assembly for new timestep. /// Preprocessing before starting assembly for new timestep.
virtual void preTimestep(GlobalVector const& /*x*/, virtual void preTimestep(GlobalVector const& /*x*/,
const double /*t*/, const double /*delta_t*/) {} const double /*t*/, const double /*delta_t*/) {}
/// Postprocessing after a complete timestep. /// Postprocessing after a complete timestep.
virtual void postTimestep(GlobalVector const& /*x*/) {} virtual void postTimestep(GlobalVector const& /*x*/) {}
/// Process output. /// Process output.
/// The file_name is indicating the name of possible output file. /// The file_name is indicating the name of possible output file.
void output(std::string const& file_name, void output(std::string const& file_name,
const unsigned /*timestep*/, const unsigned /*timestep*/,
GlobalVector const& x) const GlobalVector const& x) const
{ {
doProcessOutput(file_name, x, _mesh, *_local_to_global_index_map, doProcessOutput(file_name, x, _mesh, *_local_to_global_index_map,
_process_variables, _secondary_variables, _process_output); _process_variables, _secondary_variables, _process_output);
} }
void initialize() void initialize()
{ {
DBUG("Initialize process."); DBUG("Initialize process.");
DBUG("Construct dof mappings."); DBUG("Construct dof mappings.");
constructDofTable(); constructDofTable();
#ifndef USE_PETSC #ifndef USE_PETSC
DBUG("Compute sparsity pattern"); DBUG("Compute sparsity pattern");
computeSparsityPattern(); computeSparsityPattern();
#endif #endif
initializeConcreteProcess(*_local_to_global_index_map, _mesh, initializeConcreteProcess(*_local_to_global_index_map, _mesh,
_integration_order); _integration_order);
DBUG("Initialize boundary conditions."); DBUG("Initialize boundary conditions.");
for (ProcessVariable& pv : _process_variables) for (ProcessVariable& pv : _process_variables)
{ {
createDirichletBcs(pv, 0); // 0 is the component id createDirichletBcs(pv, 0); // 0 is the component id
createNeumannBcs(pv, 0); // 0 is the component id createNeumannBcs(pv, 0); // 0 is the component id
} }
for (auto& bc : _neumann_bcs) for (auto& bc : _neumann_bcs)
bc->initialize(_mesh.getDimension()); bc->initialize(_mesh.getDimension());
} }
void setInitialConditions(GlobalVector& x) void setInitialConditions(GlobalVector& x)
{ {
DBUG("Set initial conditions."); DBUG("Set initial conditions.");
for (ProcessVariable& pv : _process_variables) for (ProcessVariable& pv : _process_variables)
{ {
setInitialConditions(pv, 0, x); // 0 is the component id setInitialConditions(pv, 0, x); // 0 is the component id
} }
} }
MathLib::MatrixSpecifications getMatrixSpecifications() const override final MathLib::MatrixSpecifications getMatrixSpecifications() const override final
{ {
return { 0u, 0u, return { 0u, 0u,
&_sparsity_pattern, &_sparsity_pattern,
_local_to_global_index_map.get(), _local_to_global_index_map.get(),
&_mesh }; &_mesh };
} }
void assemble(const double t, GlobalVector const& x, void assemble(const double t, GlobalVector const& x,
GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override final GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override final
{ {
assembleConcreteProcess(t, x, M, K, b); assembleConcreteProcess(t, x, M, K, b);
// Call global assembler for each Neumann boundary local assembler. // Call global assembler for each Neumann boundary local assembler.
for (auto const& bc : _neumann_bcs) for (auto const& bc : _neumann_bcs)
bc->integrate(t, b); bc->integrate(t, b);
} }
void assembleJacobian( void assembleJacobian(
const double t, GlobalVector const& x, GlobalVector const& xdot, const double t, GlobalVector const& x, GlobalVector const& xdot,
const double dxdot_dx, GlobalMatrix const& M, const double dxdot_dx, GlobalMatrix const& M,
const double dx_dx, GlobalMatrix const& K, const double dx_dx, GlobalMatrix const& K,
GlobalMatrix& Jac) override final GlobalMatrix& Jac) override final
{ {
assembleJacobianConcreteProcess(t, x, xdot, dxdot_dx, M, dx_dx, K, Jac); assembleJacobianConcreteProcess(t, x, xdot, dxdot_dx, M, dx_dx, K, Jac);
// TODO In this method one could check if the user wants to use an // TODO In this method one could check if the user wants to use an
// analytical or a numerical Jacobian. Then the right // analytical or a numerical Jacobian. Then the right
// assembleJacobianConcreteProcess() method will be chosen. // assembleJacobianConcreteProcess() method will be chosen.
// Additionally in the default implementation of said method one // Additionally in the default implementation of said method one
// could provide a fallback to a numerical Jacobian. However, that // could provide a fallback to a numerical Jacobian. However, that
// would be in a sense implicit behaviour and it might be better to // would be in a sense implicit behaviour and it might be better to
// just abort, as is currently the case. // just abort, as is currently the case.
// In order to implement the Jacobian assembly entirely, in addition // In order to implement the Jacobian assembly entirely, in addition
// to the operator() in VectorMatrixAssembler there has to be a method // to the operator() in VectorMatrixAssembler there has to be a method
// that dispatches the Jacobian assembly. // that dispatches the Jacobian assembly.
// Similarly, then the NeumannBC specialization of VectorMatrixAssembler // Similarly, then the NeumannBC specialization of VectorMatrixAssembler
// probably can be merged into the main class s.t. one has only one // probably can be merged into the main class s.t. one has only one
// type of VectorMatrixAssembler (for each equation type) with the // type of VectorMatrixAssembler (for each equation type) with the
// three methods assemble(), assembleJacobian() and assembleNeumannBC(). // three methods assemble(), assembleJacobian() and assembleNeumannBC().
// That list can be extended, e.g. by methods for the assembly of // That list can be extended, e.g. by methods for the assembly of
// source terms. // source terms.
// UPDATE: Probably it is better to keep a separate NeumannBC version of the // UPDATE: Probably it is better to keep a separate NeumannBC version of the
// VectoMatrixAssembler since that will work for all kinds of processes. // VectoMatrixAssembler since that will work for all kinds of processes.
} }
std::vector<DirichletBc<Index>> const* getKnownSolutions( std::vector<DirichletBc<Index>> const* getKnownSolutions(
double const /*t*/) const override final double const /*t*/) const override final
{ {
return &_dirichlet_bcs; return &_dirichlet_bcs;
} }
NonlinearSolver& getNonlinearSolver() const NonlinearSolver& getNonlinearSolver() const
{ {
return _nonlinear_solver; return _nonlinear_solver;
} }
TimeDiscretization& getTimeDiscretization() const TimeDiscretization& getTimeDiscretization() const
{ {
return *_time_discretization; return *_time_discretization;
} }
private: private:
/// Process specific initialization called by initialize(). /// Process specific initialization called by initialize().
virtual void initializeConcreteProcess( virtual void initializeConcreteProcess(
AssemblerLib::LocalToGlobalIndexMap const& dof_table, AssemblerLib::LocalToGlobalIndexMap const& dof_table,
MeshLib::Mesh const& mesh, MeshLib::Mesh const& mesh,
unsigned const integration_order) = 0; unsigned const integration_order) = 0;
virtual void assembleConcreteProcess( virtual void assembleConcreteProcess(
const double t, GlobalVector const& x, const double t, GlobalVector const& x,
GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) = 0; GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) = 0;
virtual void assembleJacobianConcreteProcess( virtual void assembleJacobianConcreteProcess(
const double /*t*/, GlobalVector const& /*x*/, GlobalVector const& /*xdot*/, const double /*t*/, GlobalVector const& /*x*/, GlobalVector const& /*xdot*/,
const double /*dxdot_dx*/, GlobalMatrix const& /*M*/, const double /*dxdot_dx*/, GlobalMatrix const& /*M*/,
const double /*dx_dx*/, GlobalMatrix const& /*K*/, const double /*dx_dx*/, GlobalMatrix const& /*K*/,
GlobalMatrix& /*Jac*/) GlobalMatrix& /*Jac*/)
{ {
ERR("The concrete implementation of this Process did not override the" ERR("The concrete implementation of this Process did not override the"
" assembleJacobianConcreteProcess() method." " assembleJacobianConcreteProcess() method."
" Hence, no analytical Jacobian is provided for this process" " Hence, no analytical Jacobian is provided for this process"
" and the Newton-Raphson method cannot be used to solve it."); " and the Newton-Raphson method cannot be used to solve it.");
std::abort(); std::abort();
} }
void constructDofTable() void constructDofTable()
{ {
// Create single component dof in every of the mesh's nodes. // Create single component dof in every of the mesh's nodes.
_mesh_subset_all_nodes.reset( _mesh_subset_all_nodes.reset(
new MeshLib::MeshSubset(_mesh, &_mesh.getNodes())); new MeshLib::MeshSubset(_mesh, &_mesh.getNodes()));
// Collect the mesh subsets in a vector. // Collect the mesh subsets in a vector.
std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets; std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets;
for (ProcessVariable const& pv : _process_variables) for (ProcessVariable const& pv : _process_variables)
{ {
std::generate_n( std::generate_n(
std::back_inserter(all_mesh_subsets), std::back_inserter(all_mesh_subsets),
pv.getNumberOfComponents(), pv.getNumberOfComponents(),
[&]() [&]()
{ {
return std::unique_ptr<MeshLib::MeshSubsets>{ return std::unique_ptr<MeshLib::MeshSubsets>{
new MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()}}; new MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()}};
}); });
} }
_local_to_global_index_map.reset( _local_to_global_index_map.reset(
new AssemblerLib::LocalToGlobalIndexMap( new AssemblerLib::LocalToGlobalIndexMap(
std::move(all_mesh_subsets), std::move(all_mesh_subsets),
AssemblerLib::ComponentOrder::BY_LOCATION)); AssemblerLib::ComponentOrder::BY_LOCATION));
} }
/// Sets the initial condition values in the solution vector x for a given /// Sets the initial condition values in the solution vector x for a given
/// process variable and component. /// process variable and component.
void setInitialConditions(ProcessVariable const& variable, void setInitialConditions(ProcessVariable const& variable,
int const component_id, int const component_id,
GlobalVector& x) GlobalVector& x)
{ {
std::size_t const n_nodes = _mesh.getNNodes(); std::size_t const n_nodes = _mesh.getNNodes();
for (std::size_t node_id = 0; node_id < n_nodes; ++node_id) for (std::size_t node_id = 0; node_id < n_nodes; ++node_id)
{ {
MeshLib::Location const l(_mesh.getID(), MeshLib::Location const l(_mesh.getID(),
MeshLib::MeshItemType::Node, node_id); MeshLib::MeshItemType::Node, node_id);
auto global_index = std::abs( auto global_index = std::abs(
_local_to_global_index_map->getGlobalIndex(l, component_id)); _local_to_global_index_map->getGlobalIndex(l, component_id));
#ifdef USE_PETSC #ifdef USE_PETSC
// The global indices of the ghost entries of the global // The global indices of the ghost entries of the global
// matrix or the global vectors need to be set as negative values // matrix or the global vectors need to be set as negative values
// for equation assembly, however the global indices start from zero. // for equation assembly, however the global indices start from zero.
// Therefore, any ghost entry with zero index is assigned an negative // Therefore, any ghost entry with zero index is assigned an negative
// value of the vector size or the matrix dimension. // value of the vector size or the matrix dimension.
// To assign the initial value for the ghost entries, // To assign the initial value for the ghost entries,
// the negative indices of the ghost entries are restored to zero. // the negative indices of the ghost entries are restored to zero.
// checked hereby. // checked hereby.
if ( global_index == x.size() ) if ( global_index == x.size() )
global_index = 0; global_index = 0;
#endif #endif
x.set(global_index, x.set(global_index,
variable.getInitialConditionValue(*_mesh.getNode(node_id), variable.getInitialConditionValue(*_mesh.getNode(node_id),
component_id)); component_id));
} }
} }
void createDirichletBcs(ProcessVariable& variable, int const component_id) void createDirichletBcs(ProcessVariable& variable, int const component_id)
{ {
MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher = MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher =
MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher( MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(
variable.getMesh()); variable.getMesh());
variable.initializeDirichletBCs(std::back_inserter(_dirichlet_bcs), variable.initializeDirichletBCs(std::back_inserter(_dirichlet_bcs),
mesh_node_searcher, mesh_node_searcher,
*_local_to_global_index_map, *_local_to_global_index_map,
component_id); component_id);
} }
void createNeumannBcs(ProcessVariable& variable, int const component_id) void createNeumannBcs(ProcessVariable& variable, int const component_id)
{ {
// Find mesh nodes. // Find mesh nodes.
MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher = MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher =
MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher( MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(
variable.getMesh()); variable.getMesh());
MeshGeoToolsLib::BoundaryElementsSearcher mesh_element_searcher( MeshGeoToolsLib::BoundaryElementsSearcher mesh_element_searcher(
variable.getMesh(), mesh_node_searcher); variable.getMesh(), mesh_node_searcher);
// Create a neumann BC for the process variable storing them in the // Create a neumann BC for the process variable storing them in the
// _neumann_bcs vector. // _neumann_bcs vector.
variable.createNeumannBcs<GlobalSetup>( variable.createNeumannBcs<GlobalSetup>(
std::back_inserter(_neumann_bcs), std::back_inserter(_neumann_bcs),
mesh_element_searcher, mesh_element_searcher,
_integration_order, _integration_order,
*_local_to_global_index_map, *_local_to_global_index_map,
0, // 0 is the variable id TODO 0, // 0 is the variable id TODO
component_id); component_id);
} }
/// Computes and stores global matrix' sparsity pattern from given /// Computes and stores global matrix' sparsity pattern from given
/// DOF-table. /// DOF-table.
void computeSparsityPattern() void computeSparsityPattern()
{ {
_sparsity_pattern = std::move(AssemblerLib::computeSparsityPattern( _sparsity_pattern = std::move(AssemblerLib::computeSparsityPattern(
*_local_to_global_index_map, _mesh)); *_local_to_global_index_map, _mesh));
} }
private: private:
unsigned const _integration_order = 2; unsigned const _integration_order = 2;
MeshLib::Mesh& _mesh; MeshLib::Mesh& _mesh;
std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap> std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap>
_local_to_global_index_map; _local_to_global_index_map;
MathLib::SparsityPattern _sparsity_pattern; MathLib::SparsityPattern _sparsity_pattern;
std::vector<DirichletBc<GlobalIndexType>> _dirichlet_bcs; std::vector<DirichletBc<GlobalIndexType>> _dirichlet_bcs;
std::vector<std::unique_ptr<NeumannBc<GlobalSetup>>> _neumann_bcs; std::vector<std::unique_ptr<NeumannBc<GlobalSetup>>> _neumann_bcs;
NonlinearSolver& _nonlinear_solver; NonlinearSolver& _nonlinear_solver;
std::unique_ptr<TimeDiscretization> _time_discretization; std::unique_ptr<TimeDiscretization> _time_discretization;
/// Variables used by this process. /// Variables used by this process.
std::vector<std::reference_wrapper<ProcessVariable>> _process_variables; std::vector<std::reference_wrapper<ProcessVariable>> _process_variables;
ProcessOutput<GlobalVector> _process_output; ProcessOutput<GlobalVector> _process_output;
protected: protected:
std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_all_nodes; std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_all_nodes;
SecondaryVariableCollection<GlobalVector> _secondary_variables; SecondaryVariableCollection<GlobalVector> _secondary_variables;
}; };
/// Find process variables in \c variables whose names match the settings under /// Find process variables in \c variables whose names match the settings under
...@@ -366,36 +366,36 @@ Parameter<ParameterArgs...>& findParameter( ...@@ -366,36 +366,36 @@ Parameter<ParameterArgs...>& findParameter(
BaseLib::ConfigTree const& process_config, std::string const& tag, BaseLib::ConfigTree const& process_config, std::string const& tag,
std::vector<std::unique_ptr<ParameterBase>> const& parameters) std::vector<std::unique_ptr<ParameterBase>> const& parameters)
{ {
// Find parameter name in process config. // Find parameter name in process config.
auto const name = process_config.getConfParam<std::string>(tag); auto const name = process_config.getConfParam<std::string>(tag);
// Find corresponding parameter by name. // Find corresponding parameter by name.
auto const parameter_it = auto const parameter_it =
std::find_if(parameters.cbegin(), parameters.cend(), std::find_if(parameters.cbegin(), parameters.cend(),
[&name](std::unique_ptr<ParameterBase> const& p) [&name](std::unique_ptr<ParameterBase> const& p)
{ {
return p->name == name; return p->name == name;
}); });
if (parameter_it == parameters.end()) if (parameter_it == parameters.end())
{ {
ERR( ERR(
"Could not find parameter '%s' in the provided parameters list for " "Could not find parameter '%s' in the provided parameters list for "
"config tag <%s>.", "config tag <%s>.",
name.c_str(), tag.c_str()); name.c_str(), tag.c_str());
std::abort(); std::abort();
} }
DBUG("Found parameter \'%s\'.", (*parameter_it)->name.c_str()); DBUG("Found parameter \'%s\'.", (*parameter_it)->name.c_str());
// Check the type correctness of the found parameter. // Check the type correctness of the found parameter.
auto* const parameter = auto* const parameter =
dynamic_cast<Parameter<ParameterArgs...>*>(parameter_it->get()); dynamic_cast<Parameter<ParameterArgs...>*>(parameter_it->get());
if (!parameter) if (!parameter)
{ {
ERR("The read parameter is of incompatible type."); ERR("The read parameter is of incompatible type.");
std::abort(); std::abort();
} }
return *parameter; return *parameter;
} }
} // namespace ProcessLib } // namespace ProcessLib
......
...@@ -25,74 +25,74 @@ ProcessVariable::ProcessVariable(BaseLib::ConfigTree const& config, ...@@ -25,74 +25,74 @@ ProcessVariable::ProcessVariable(BaseLib::ConfigTree const& config,
_mesh(mesh), _mesh(mesh),
_n_components(config.getConfParam<int>("components")) _n_components(config.getConfParam<int>("components"))
{ {
DBUG("Constructing process variable %s", this->_name.c_str()); DBUG("Constructing process variable %s", this->_name.c_str());
// Initial condition // Initial condition
if (auto ic_config = config.getConfSubtreeOptional("initial_condition")) if (auto ic_config = config.getConfSubtreeOptional("initial_condition"))
{ {
auto const type = ic_config->peekConfParam<std::string>("type"); auto const type = ic_config->peekConfParam<std::string>("type");
if (type == "Uniform") if (type == "Uniform")
{ {
_initial_condition = _initial_condition =
createUniformInitialCondition(*ic_config, _n_components); createUniformInitialCondition(*ic_config, _n_components);
} }
else if (type == "MeshProperty") else if (type == "MeshProperty")
{ {
_initial_condition = _initial_condition =
createMeshPropertyInitialCondition(*ic_config, _mesh, _n_components); createMeshPropertyInitialCondition(*ic_config, _mesh, _n_components);
} }
else else
{ {
ERR("Unknown type of the initial condition."); ERR("Unknown type of the initial condition.");
} }
} }
else else
{ {
INFO("No initial condition found."); INFO("No initial condition found.");
} }
// Boundary conditions // Boundary conditions
if (auto bcs_config = config.getConfSubtreeOptional("boundary_conditions")) if (auto bcs_config = config.getConfSubtreeOptional("boundary_conditions"))
{ {
for (auto bc_config for (auto bc_config
: bcs_config->getConfSubtreeList("boundary_condition")) : bcs_config->getConfSubtreeList("boundary_condition"))
{ {
auto const geometrical_set_name = auto const geometrical_set_name =
bc_config.getConfParam<std::string>("geometrical_set"); bc_config.getConfParam<std::string>("geometrical_set");
auto const geometry_name = auto const geometry_name =
bc_config.getConfParam<std::string>("geometry"); bc_config.getConfParam<std::string>("geometry");
GeoLib::GeoObject const* const geometry = GeoLib::GeoObject const* const geometry =
geometries.getGeoObject(geometrical_set_name, geometry_name); geometries.getGeoObject(geometrical_set_name, geometry_name);
DBUG( DBUG(
"Found geometry type \"%s\"", "Found geometry type \"%s\"",
GeoLib::convertGeoTypeToString(geometry->getGeoType()).c_str()); GeoLib::convertGeoTypeToString(geometry->getGeoType()).c_str());
// Construct type dependent boundary condition // Construct type dependent boundary condition
auto const type = bc_config.peekConfParam<std::string>("type"); auto const type = bc_config.peekConfParam<std::string>("type");
if (type == "UniformDirichlet") if (type == "UniformDirichlet")
{ {
_dirichlet_bc_configs.emplace_back( _dirichlet_bc_configs.emplace_back(
new UniformDirichletBoundaryCondition(geometry, bc_config)); new UniformDirichletBoundaryCondition(geometry, bc_config));
} }
else if (type == "UniformNeumann") else if (type == "UniformNeumann")
{ {
_neumann_bc_configs.emplace_back( _neumann_bc_configs.emplace_back(
new NeumannBcConfig(geometry, bc_config)); new NeumannBcConfig(geometry, bc_config));
} }
else else
{ {
ERR("Unknown type \'%s\' of the boundary condition.", ERR("Unknown type \'%s\' of the boundary condition.",
type.c_str()); type.c_str());
} }
} }
} else { } else {
INFO("No boundary conditions found."); INFO("No boundary conditions found.");
} }
// Source Terms // Source Terms
config.ignoreConfParam("source_terms"); config.ignoreConfParam("source_terms");
} }
ProcessVariable::ProcessVariable(ProcessVariable&& other) ProcessVariable::ProcessVariable(ProcessVariable&& other)
...@@ -107,32 +107,32 @@ ProcessVariable::ProcessVariable(ProcessVariable&& other) ...@@ -107,32 +107,32 @@ ProcessVariable::ProcessVariable(ProcessVariable&& other)
std::string const& ProcessVariable::getName() const std::string const& ProcessVariable::getName() const
{ {
return _name; return _name;
} }
MeshLib::Mesh const& ProcessVariable::getMesh() const MeshLib::Mesh const& ProcessVariable::getMesh() const
{ {
return _mesh; return _mesh;
} }
MeshLib::PropertyVector<double>& ProcessVariable::getOrCreateMeshProperty() MeshLib::PropertyVector<double>& ProcessVariable::getOrCreateMeshProperty()
{ {
boost::optional<MeshLib::PropertyVector<double>&> result; boost::optional<MeshLib::PropertyVector<double>&> result;
if (_mesh.getProperties().hasPropertyVector(_name)) if (_mesh.getProperties().hasPropertyVector(_name))
{ {
result = result =
_mesh.getProperties().template getPropertyVector<double>(_name); _mesh.getProperties().template getPropertyVector<double>(_name);
assert(result); assert(result);
assert(result->size() == _mesh.getNNodes() * _n_components); assert(result->size() == _mesh.getNNodes() * _n_components);
} }
else else
{ {
result = _mesh.getProperties().template createNewPropertyVector<double>( result = _mesh.getProperties().template createNewPropertyVector<double>(
_name, MeshLib::MeshItemType::Node); _name, MeshLib::MeshItemType::Node);
assert(result); assert(result);
result->resize(_mesh.getNNodes() * _n_components); result->resize(_mesh.getNNodes() * _n_components);
} }
return *result; return *result;
} }
} // namespace ProcessLib } // namespace ProcessLib
...@@ -47,67 +47,67 @@ namespace ProcessLib ...@@ -47,67 +47,67 @@ namespace ProcessLib
class ProcessVariable class ProcessVariable
{ {
public: public:
ProcessVariable(BaseLib::ConfigTree const& config, MeshLib::Mesh& mesh, ProcessVariable(BaseLib::ConfigTree const& config, MeshLib::Mesh& mesh,
GeoLib::GEOObjects const& geometries); GeoLib::GEOObjects const& geometries);
ProcessVariable(ProcessVariable&&); ProcessVariable(ProcessVariable&&);
std::string const& getName() const; std::string const& getName() const;
/// Returns a mesh on which the process variable is defined. /// Returns a mesh on which the process variable is defined.
MeshLib::Mesh const& getMesh() const; MeshLib::Mesh const& getMesh() const;
/// Returns the number of components of the process variable. /// Returns the number of components of the process variable.
int getNumberOfComponents() const { return _n_components; } int getNumberOfComponents() const { return _n_components; }
template <typename OutputIterator> template <typename OutputIterator>
void initializeDirichletBCs( void initializeDirichletBCs(
OutputIterator output_bcs, OutputIterator output_bcs,
MeshGeoToolsLib::MeshNodeSearcher& searcher, MeshGeoToolsLib::MeshNodeSearcher& searcher,
const AssemblerLib::LocalToGlobalIndexMap& dof_table, const AssemblerLib::LocalToGlobalIndexMap& dof_table,
const unsigned nodal_dof_idx) const unsigned nodal_dof_idx)
{ {
for (auto& bc_config : _dirichlet_bc_configs) for (auto& bc_config : _dirichlet_bc_configs)
{ {
DirichletBc<GlobalIndexType> bc; DirichletBc<GlobalIndexType> bc;
bc_config->initialize(searcher, dof_table, nodal_dof_idx, bc); bc_config->initialize(searcher, dof_table, nodal_dof_idx, bc);
output_bcs++ = bc; output_bcs++ = bc;
} }
} }
template <typename GlobalSetup, typename OutputIterator, typename... Args> template <typename GlobalSetup, typename OutputIterator, typename... Args>
void createNeumannBcs(OutputIterator bcs, void createNeumannBcs(OutputIterator bcs,
MeshGeoToolsLib::BoundaryElementsSearcher& searcher, MeshGeoToolsLib::BoundaryElementsSearcher& searcher,
Args&&... args) Args&&... args)
{ {
for (auto& config : _neumann_bc_configs) for (auto& config : _neumann_bc_configs)
{ {
config->initialize(searcher); config->initialize(searcher);
bcs++ = std::unique_ptr<NeumannBc<GlobalSetup>>{ bcs++ = std::unique_ptr<NeumannBc<GlobalSetup>>{
new NeumannBc<GlobalSetup>(*config, new NeumannBc<GlobalSetup>(*config,
std::forward<Args>(args)...)}; std::forward<Args>(args)...)};
} }
} }
double getInitialConditionValue(MeshLib::Node const& n, double getInitialConditionValue(MeshLib::Node const& n,
int const component_id) const int const component_id) const
{ {
return _initial_condition->getValue(n, component_id); return _initial_condition->getValue(n, component_id);
} }
// Get or create a property vector for results. // Get or create a property vector for results.
// The returned mesh property size is number of mesh nodes times number of // The returned mesh property size is number of mesh nodes times number of
// components. // components.
MeshLib::PropertyVector<double>& getOrCreateMeshProperty(); MeshLib::PropertyVector<double>& getOrCreateMeshProperty();
private: private:
std::string const _name; std::string const _name;
MeshLib::Mesh& _mesh; MeshLib::Mesh& _mesh;
const int _n_components; const int _n_components;
std::unique_ptr<InitialCondition> _initial_condition; std::unique_ptr<InitialCondition> _initial_condition;
std::vector<std::unique_ptr<UniformDirichletBoundaryCondition>> std::vector<std::unique_ptr<UniformDirichletBoundaryCondition>>
_dirichlet_bc_configs; _dirichlet_bc_configs;
std::vector<std::unique_ptr<NeumannBcConfig>> _neumann_bc_configs; std::vector<std::unique_ptr<NeumannBcConfig>> _neumann_bc_configs;
}; };
} // namespace ProcessLib } // namespace ProcessLib
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment