Skip to content

Commit

Permalink
Support a different elastic material in each block
Browse files Browse the repository at this point in the history
  • Loading branch information
nilsvu committed Feb 29, 2024
1 parent 5bd55a0 commit 670e39c
Show file tree
Hide file tree
Showing 20 changed files with 566 additions and 78 deletions.
1 change: 1 addition & 0 deletions src/Elliptic/Executables/Elasticity/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ set(LIBS_TO_LINK
DiscontinuousGalerkin
DomainCreators
Elasticity
ElasticityActions
ElasticityBoundaryConditions
ElasticityPointwiseFunctions
ElasticitySolutions
Expand Down
25 changes: 17 additions & 8 deletions src/Elliptic/Executables/Elasticity/SolveElasticity.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include "Elliptic/BoundaryConditions/BoundaryCondition.hpp"
#include "Elliptic/DiscontinuousGalerkin/DgElementArray.hpp"
#include "Elliptic/Executables/Solver.hpp"
#include "Elliptic/Systems/Elasticity/Actions/InitializeConstitutiveRelation.hpp"
#include "Elliptic/Systems/Elasticity/BoundaryConditions/Factory.hpp"
#include "Elliptic/Systems/Elasticity/FirstOrderSystem.hpp"
#include "Elliptic/Systems/Elasticity/Tags.hpp"
Expand Down Expand Up @@ -78,8 +79,10 @@ struct Metavariables {
error_compute>;

// Collect all items to store in the cache.
using const_global_cache_tags =
tmpl::list<Elasticity::Tags::ConstitutiveRelation<volume_dim>>;
using const_global_cache_tags = tmpl::list<>;

// Just a name in the input file
struct ObserveNormsPerLayer {};

struct factory_creation
: tt::ConformsTo<Options::protocols::FactoryCreation> {
Expand All @@ -98,12 +101,17 @@ struct Metavariables {
Elasticity::ConstitutiveRelations::ConstitutiveRelation<volume_dim>,
Elasticity::ConstitutiveRelations::standard_constitutive_relations<
volume_dim>>,
tmpl::pair<Event,
tmpl::flatten<tmpl::list<
Events::Completion,
dg::Events::field_observations<
volume_dim, observe_fields, observer_compute_tags,
LinearSolver::multigrid::Tags::IsFinestGrid>>>>,
tmpl::pair<
Event,
tmpl::flatten<tmpl::list<
Events::Completion,
dg::Events::field_observations<
volume_dim, observe_fields, observer_compute_tags,
LinearSolver::multigrid::Tags::IsFinestGrid>,
// Observation per material layer
::Events::ObserveNorms<observe_fields, observer_compute_tags,
Elasticity::Tags::MaterialLayerName,
ObserveNormsPerLayer>>>>,
tmpl::pair<Trigger, elliptic::Triggers::all_triggers<
typename solver::linear_solver::options_group>>,
tmpl::pair<PhaseChange, tmpl::list<PhaseControl::VisitAndReturn<
Expand All @@ -118,6 +126,7 @@ struct Metavariables {

using initialization_actions =
tmpl::push_back<typename solver::initialization_actions,
Elasticity::Actions::InitializeConstitutiveRelation<Dim>,
Parallel::Actions::TerminatePhase>;

using register_actions =
Expand Down
31 changes: 31 additions & 0 deletions src/Elliptic/Systems/Elasticity/Actions/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
# Distributed under the MIT License.
# See LICENSE.txt for details.

set(LIBRARY ElasticityActions)

add_spectre_library(${LIBRARY})

spectre_target_sources(
${LIBRARY}
PRIVATE
InitializeConstitutiveRelation.cpp
)

spectre_target_headers(
${LIBRARY}
INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
HEADERS
InitializeConstitutiveRelation.hpp
)

target_link_libraries(
${LIBRARY}
PUBLIC
ConstitutiveRelations
DomainCreators
INTERFACE
DataStructures
Domain
Parallel
Utilities
)
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
// Distributed under the MIT License.
// See LICENSE.txt for details.

// Instantiations of domain::ExpandOverBlocks for Elasticity

#include <cstddef>
#include <memory>

#include "Domain/Creators/ExpandOverBlocks.tpp"
#include "PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp"

template <size_t Dim>
using ConstRelPtr = std::unique_ptr<
Elasticity::ConstitutiveRelations::ConstitutiveRelation<Dim>>;

template class domain::ExpandOverBlocks<ConstRelPtr<2>>;
template class domain::ExpandOverBlocks<ConstRelPtr<3>>;
Original file line number Diff line number Diff line change
@@ -0,0 +1,201 @@
// Distributed under the MIT License.
// See LICENSE.txt for details.

#pragma once

#include <cstddef>
#include <exception>
#include <memory>
#include <string>
#include <tuple>
#include <unordered_map>
#include <utility>
#include <variant>
#include <vector>

#include "DataStructures/DataBox/DataBox.hpp"
#include "DataStructures/DataBox/Tag.hpp"
#include "Domain/Creators/BlockGroups.hpp"
#include "Domain/Creators/DomainCreator.hpp"
#include "Domain/Creators/ExpandOverBlocks.hpp"
#include "Domain/Creators/OptionTags.hpp"
#include "Domain/Creators/Tags/Domain.hpp"
#include "Domain/Structure/ElementId.hpp"
#include "Domain/Tags.hpp"
#include "Elliptic/Tags.hpp"
#include "IO/Observer/Tags.hpp"
#include "Options/String.hpp"
#include "Parallel/AlgorithmExecution.hpp"
#include "PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp"
#include "PointwiseFunctions/Elasticity/ConstitutiveRelations/Tags.hpp"
#include "Utilities/CallWithDynamicType.hpp"
#include "Utilities/ErrorHandling/Error.hpp"
#include "Utilities/TMPL.hpp"

/// \cond
namespace tuples {
template <typename... Tags>
struct TaggedTuple;
} // namespace tuples
namespace Parallel {
template <typename Metavariables>
struct GlobalCache;
} // namespace Parallel
/// \endcond

namespace Elasticity {
namespace Tags {

/// A constitutive relation in every block of the domain
template <size_t Dim>
struct ConstitutiveRelationPerBlock : db::SimpleTag,
ConstitutiveRelationPerBlockBase {
using ConstRelPtr =
std::unique_ptr<ConstitutiveRelations::ConstitutiveRelation<Dim>>;
using type = std::vector<ConstRelPtr>;

using option_tags = tmpl::list<domain::OptionTags::DomainCreator<Dim>,
OptionTags::ConstitutiveRelationPerBlock<Dim>>;
static constexpr bool pass_metavariables = false;

static type create_from_options(
const std::unique_ptr<DomainCreator<Dim>>& domain_creator,
const std::variant<ConstRelPtr, std::vector<ConstRelPtr>,
std::unordered_map<std::string, ConstRelPtr>>&
constitutive_relation_per_block) {
const auto block_names = domain_creator->block_names();
const auto block_groups = domain_creator->block_groups();
const domain::ExpandOverBlocks<ConstRelPtr> expand_over_blocks{
block_names, block_groups};
try {
return std::visit(expand_over_blocks, constitutive_relation_per_block);
} catch (const std::exception& error) {
ERROR_NO_TRACE("Invalid 'Material':\n" << error.what());
}
}
};

/// References the constitutive relation for the element's block, which is
/// stored in the global cache
template <size_t Dim>
struct ConstitutiveRelationReference : ConstitutiveRelation<Dim>,
db::ReferenceTag {
using base = ConstitutiveRelation<Dim>;
using argument_tags =
tmpl::list<ConstitutiveRelationPerBlockBase, domain::Tags::Element<Dim>>;
static const ConstitutiveRelations::ConstitutiveRelation<Dim>& get(
const std::vector<
std::unique_ptr<ConstitutiveRelations::ConstitutiveRelation<Dim>>>&
constitutive_relation_per_block,
const Element<Dim>& element) {
return *constitutive_relation_per_block.at(element.id().block_id());
}
};

/// Stores the names of the block groups that split the domain into layers with
/// different material properties. Useful to observe quantities in each layer.
template <size_t Dim>
struct MaterialBlockGroups : db::SimpleTag {
using type = std::unordered_set<std::string>;

using option_tags = tmpl::list<OptionTags::ConstitutiveRelationPerBlock<Dim>>;
static constexpr bool pass_metavariables = false;
using ConstRelPtr =
std::unique_ptr<ConstitutiveRelations::ConstitutiveRelation<Dim>>;

static type create_from_options(
const std::variant<ConstRelPtr, std::vector<ConstRelPtr>,
std::unordered_map<std::string, ConstRelPtr>>&
constitutive_relation_per_block) {
if (std::holds_alternative<std::unordered_map<std::string, ConstRelPtr>>(
constitutive_relation_per_block)) {
const auto& map = std::get<std::unordered_map<std::string, ConstRelPtr>>(
constitutive_relation_per_block);
std::unordered_set<std::string> block_groups;
for (const auto& [block_name, _] : map) {
block_groups.insert(block_name);
}
return block_groups;
} else {
return {};
}
}
};

/// The name of the material layer (name of a block group with some material)
struct MaterialLayerName : db::SimpleTag {
using type = std::optional<std::string>;
};

} // namespace Tags

/// Actions related to solving Elasticity systems
namespace Actions {

/*!
* \brief Initialize the constitutive relation describing properties of the
* elastic material
*
* Every block in the domain can have a different constitutive relation,
* allowing for composite materials. All constitutive relations are stored in
* the global cache indexed by block, and elements reference their block's
* constitutive relation in the DataBox. This means an element can retrieve the
* local constitutive relation from the DataBox simply by requesting
* `Elasticity::Tags::ConstitutiveRelation<Dim>`.
*/
template <size_t Dim>
struct InitializeConstitutiveRelation {
public:
using const_global_cache_tags =
tmpl::list<Tags::ConstitutiveRelationPerBlock<Dim>,
Tags::MaterialBlockGroups<Dim>>;
using simple_tags =
tmpl::list<Tags::MaterialLayerName,
observers::Tags::ObservationKey<Tags::MaterialLayerName>>;
using compute_tags = tmpl::list<Tags::ConstitutiveRelationReference<Dim>>;

template <typename DbTags, typename... InboxTags, typename Metavariables,
typename ActionList, typename ParallelComponent>
static Parallel::iterable_action_return_t apply(
db::DataBox<DbTags>& box,
const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
const Parallel::GlobalCache<Metavariables>& /*cache*/,
const ElementId<Dim>& /*array_index*/, const ActionList /*meta*/,
const ParallelComponent* const /*meta*/) {
db::mutate_apply<InitializeConstitutiveRelation>(make_not_null(&box));
return {Parallel::AlgorithmExecution::Continue, std::nullopt};
}

public:
using return_tags = simple_tags;
using argument_tags =
tmpl::list<Tags::MaterialBlockGroups<Dim>, domain::Tags::Element<Dim>,
domain::Tags::Domain<Dim>>;

static void apply(
const gsl::not_null<std::optional<std::string>*> material_layer_name,
const gsl::not_null<std::optional<std::string>*> observation_key,
const std::unordered_set<std::string>& material_block_groups,
const Element<Dim>& element, const Domain<Dim>& domain) {
const auto& block = domain.blocks()[element.id().block_id()];
// Check if this element is in a material layer
*material_layer_name = [&material_block_groups, &domain,
&block]() -> std::optional<std::string> {
for (const auto& name : material_block_groups) {
if (domain::block_is_in_group(block.name(), name,
domain.block_groups())) {
return name;
}
}
return std::nullopt;
}();
// Set the corresponding observation key, but only on the finest multigrid
// level. This could be done better by supporting intersections of array
// sections in observation events or something like that.
*observation_key =
element.id().grid_index() == 0 ? *material_layer_name : std::nullopt;
}
};

} // namespace Actions
} // namespace Elasticity
1 change: 1 addition & 0 deletions src/Elliptic/Systems/Elasticity/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -31,4 +31,5 @@ target_link_libraries(
LinearOperators
)

add_subdirectory(Actions)
add_subdirectory(BoundaryConditions)
21 changes: 14 additions & 7 deletions src/Elliptic/Systems/Elasticity/Equations.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@

#include <algorithm>
#include <cstddef>
#include <memory>
#include <vector>

#include "DataStructures/DataBox/PrefixHelpers.hpp"
#include "DataStructures/DataBox/Prefixes.hpp"
Expand Down Expand Up @@ -57,21 +59,24 @@ void add_curved_sources(
template <size_t Dim>
void Fluxes<Dim>::apply(
const gsl::not_null<tnsr::II<DataVector, Dim>*> minus_stress,
const ConstitutiveRelations::ConstitutiveRelation<Dim>&
constitutive_relation,
const tnsr::I<DataVector, Dim>& coordinates,
const std::vector<
std::unique_ptr<ConstitutiveRelations::ConstitutiveRelation<Dim>>>&
constitutive_relation_per_block,
const Element<Dim>& element, const tnsr::I<DataVector, Dim>& coordinates,
const tnsr::I<DataVector, Dim>& /*displacement*/,
const tnsr::iJ<DataVector, Dim>& deriv_displacement) {
primal_fluxes(minus_stress, deriv_displacement, constitutive_relation,
primal_fluxes(minus_stress, deriv_displacement,
*constitutive_relation_per_block.at(element.id().block_id()),
coordinates);
}

template <size_t Dim>
void Fluxes<Dim>::apply(
const gsl::not_null<tnsr::II<DataVector, Dim>*> minus_stress,
const ConstitutiveRelations::ConstitutiveRelation<Dim>&
constitutive_relation,
const tnsr::I<DataVector, Dim>& coordinates,
const std::vector<
std::unique_ptr<ConstitutiveRelations::ConstitutiveRelation<Dim>>>&
constitutive_relation_per_block,
const Element<Dim>& element, const tnsr::I<DataVector, Dim>& coordinates,
const tnsr::i<DataVector, Dim>& face_normal,
const tnsr::I<DataVector, Dim>& /*face_normal_vector*/,
const tnsr::I<DataVector, Dim>& displacement) {
Expand All @@ -82,6 +87,8 @@ void Fluxes<Dim>::apply(
face_normal.get(j) * displacement.get(i));
}
}
const auto& constitutive_relation =
*constitutive_relation_per_block.at(element.id().block_id());
constitutive_relation.stress(minus_stress, strain, coordinates);
for (auto& component : *minus_stress) {
component *= -1.;
Expand Down
Loading

0 comments on commit 670e39c

Please sign in to comment.