Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add domain of definition for parameters. #2376

Merged
merged 12 commits into from Mar 1, 2019

[PL] BC; Cleanup assemblies after switch to params

  • Loading branch information...
endJunction committed Feb 22, 2019
commit c62b83dd014ed2b4d4123284d51372ec5081454d
@@ -54,28 +54,14 @@ class NeumannBoundaryConditionLocalAssembler final
unsigned const n_integration_points =
Base::_integration_method.getNumberOfPoints();

SpatialPosition pos;
pos.setElementID(id);

using FemType =
NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
FemType fe(
static_cast<const typename ShapeFunction::MeshElement&>(Base::_element));

// Get element nodes for the interpolation from nodes to
// integration point.
// Get element nodes for the interpolation from nodes to integration
// point.
NodalVectorType parameter_node_values =
_neumann_bc_parameter.getNodalValuesOnElement(Base::_element, t);

for (unsigned ip = 0; ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);
auto const& ip_data = Base::_ns_and_weights[ip];
auto const& wp = Base::_integration_method.getWeightedPoint(ip);

MathLib::TemplatePoint<double, 3> coords_ip(
fe.interpolateCoordinates(ip_data.N));
pos.setCoordinates(coords_ip);

_local_rhs.noalias() += ip_data.N *
parameter_node_values.dot(ip_data.N) *
@@ -21,21 +21,19 @@ namespace ProcessLib
{
namespace NormalTractionBoundaryCondition
{
template <typename ShapeMatricesTypeDisplacement, int GlobalDim, int NPoints>
template <typename ShapeMatricesType>
struct IntegrationPointData final
{
IntegrationPointData(
typename ShapeMatricesTypeDisplacement::template VectorType<
NPoints * GlobalDim> const& Nu_times_n_,
typename ShapeMatricesType::ShapeMatrices::ShapeType const N_,
typename ShapeMatricesType::GlobalDimVectorType const n_,
double const integration_weight_)
: Nu_times_n(Nu_times_n_), integration_weight(integration_weight_)
: N(N_), n(n_), integration_weight(integration_weight_)
{
}

// Shape matrix (for displacement) times element's normal vector.
typename ShapeMatricesTypeDisplacement::template VectorType<
NPoints * GlobalDim> const Nu_times_n;

typename ShapeMatricesType::ShapeMatrices::ShapeType const N;
typename ShapeMatricesType::GlobalDimVectorType const n;
double const integration_weight;
};

@@ -56,19 +54,20 @@ class NormalTractionBoundaryConditionLocalAssembler final
: public NormalTractionBoundaryConditionLocalAssemblerInterface
{
public:
using ShapeMatricesTypeDisplacement =
using ShapeMatricesType =
ShapeMatrixPolicyType<ShapeFunctionDisplacement, GlobalDim>;
using GlobalDimVectorType =
typename ShapeMatrixPolicyType<ShapeFunctionDisplacement,
GlobalDim>::GlobalDimVectorType;
using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
using NodalVectorType = typename ShapeMatricesType::NodalVectorType;

NormalTractionBoundaryConditionLocalAssembler(
MeshLib::Element const& e,
std::size_t const local_matrix_size,
bool const is_axially_symmetric,
unsigned const integration_order,
Parameter<double> const& pressure)
: _integration_method(integration_order), _pressure(pressure)
: _integration_method(integration_order),
_pressure(pressure),
_element(e)
{
_local_rhs.setZero(local_matrix_size);

@@ -78,10 +77,9 @@ class NormalTractionBoundaryConditionLocalAssembler final
_ip_data.reserve(n_integration_points);

auto const shape_matrices_u =
initShapeMatrices<ShapeFunctionDisplacement,
ShapeMatricesTypeDisplacement, IntegrationMethod,
GlobalDim>(e, is_axially_symmetric,
_integration_method);
initShapeMatrices<ShapeFunctionDisplacement, ShapeMatricesType,
IntegrationMethod, GlobalDim>(
e, is_axially_symmetric, _integration_method);

GlobalDimVectorType element_normal(GlobalDim);

@@ -104,24 +102,13 @@ class NormalTractionBoundaryConditionLocalAssembler final

for (unsigned ip = 0; ip < n_integration_points; ip++)
{
typename ShapeMatricesTypeDisplacement::template MatrixType<
GlobalDim, displacement_size>
N_u = ShapeMatricesTypeDisplacement::template MatrixType<
GlobalDim, displacement_size>::Zero(GlobalDim,
displacement_size);
for (int i = 0; i < static_cast<int>(GlobalDim); ++i)
{
N_u.template block<1, displacement_size / GlobalDim>(
i, i * displacement_size / GlobalDim)
.noalias() = shape_matrices_u[ip].N;
}

double const integration_weight =
_integration_method.getWeightedPoint(ip).getWeight() *
shape_matrices_u[ip].integralMeasure *
shape_matrices_u[ip].detJ;

_ip_data.emplace_back(N_u.transpose() * element_normal,
_ip_data.emplace_back(shape_matrices_u[ip].N, element_normal,
integration_weight);
}
}
@@ -137,18 +124,28 @@ class NormalTractionBoundaryConditionLocalAssembler final
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();

SpatialPosition pos;
pos.setElementID(id);
NodalVectorType pressure =
_pressure.getNodalValuesOnElement(_element, t);

for (unsigned ip = 0; ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);

auto const& w = _ip_data[ip].integration_weight;
auto const& Nu_times_n = _ip_data[ip].Nu_times_n;
auto const& N = _ip_data[ip].N;
auto const& n = _ip_data[ip].n;

_local_rhs.noalias() -=
Nu_times_n.transpose() * _pressure(t, pos)[0] * w;
typename ShapeMatricesType::template MatrixType<GlobalDim,
displacement_size>
N_u = ShapeMatricesType::template MatrixType<
GlobalDim, displacement_size>::Zero(GlobalDim,
displacement_size);
for (int i = 0; i < static_cast<int>(GlobalDim); ++i)
{
N_u.template block<1, displacement_size / GlobalDim>(
i, i * displacement_size / GlobalDim)
.noalias() = N;

This comment has been minimized.

Copy link
@TomFischer

TomFischer Feb 25, 2019

Member

Curly brackets for the for-loop-body are missing. ✔️

}

_local_rhs.noalias() -= n.transpose() * N_u * pressure.dot(N) * w;
}

auto const indices = NumLib::getIndices(id, dof_table_boundary);
@@ -161,17 +158,16 @@ class NormalTractionBoundaryConditionLocalAssembler final

static const int displacement_size =
ShapeFunctionDisplacement::NPOINTS * GlobalDim;
std::vector<IntegrationPointData<ShapeMatricesTypeDisplacement, GlobalDim,
ShapeFunctionDisplacement::NPOINTS>,
Eigen::aligned_allocator<IntegrationPointData<
ShapeMatricesTypeDisplacement, GlobalDim,
ShapeFunctionDisplacement::NPOINTS>>>
std::vector<
IntegrationPointData<ShapeMatricesType>,
Eigen::aligned_allocator<IntegrationPointData<ShapeMatricesType>>>
_ip_data;

typename ShapeMatricesTypeDisplacement::template VectorType<
displacement_size>
typename ShapeMatricesType::template VectorType<displacement_size>
_local_rhs;

MeshLib::Element const& _element;

public:
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
@@ -41,9 +41,7 @@ class PythonBoundaryConditionLocalAssembler final
bool is_axially_symmetric,
unsigned const integration_order,
PythonBoundaryConditionData const& data)
: Base(e, is_axially_symmetric, integration_order),
_data(data),
_element(e)
: Base(e, is_axially_symmetric, integration_order), _data(data)
{
}

@@ -57,12 +55,12 @@ class PythonBoundaryConditionLocalAssembler final
using FemType =
NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
&_element));
&this->_element));

unsigned const num_integration_points =
Base::_integration_method.getNumberOfPoints();
auto const num_var = _data.dof_table_bulk.getNumberOfVariables();
auto const num_nodes = _element.getNumberOfNodes();
auto const num_nodes = Base::_element.getNumberOfNodes();
auto const num_comp_total =
_data.dof_table_bulk.getNumberOfComponents();

@@ -85,7 +83,8 @@ class PythonBoundaryConditionLocalAssembler final
for (unsigned element_node_id = 0; element_node_id < num_nodes;
++element_node_id)
{
auto const* const node = _element.getNode(element_node_id);
auto const* const node =
Base::_element.getNode(element_node_id);
auto const boundary_node_id = node->getID();
auto const bulk_node_id =
bulk_node_ids_map[boundary_node_id];
@@ -119,12 +118,12 @@ class PythonBoundaryConditionLocalAssembler final

for (unsigned ip = 0; ip < num_integration_points; ip++)
{
auto const& sm = Base::_shape_matrices[ip];
auto const coords = fe.interpolateCoordinates(sm.N);
auto const& N = Base::_ns_and_weights[ip].N;
auto const& w = Base::_ns_and_weights[ip].weight;
auto const coords = fe.interpolateCoordinates(N);
prim_vars =
sm.N *
primary_variables_mat; // Assumption: all primary variables
// have same shape functions.
N * primary_variables_mat; // Assumption: all primary variables
// have same shape functions.
auto const flag_flux_dFlux =
_data.bc_object->getFlux(t, coords, prim_vars_data);
if (!_data.bc_object->isOverriddenNatural())
@@ -143,9 +142,7 @@ class PythonBoundaryConditionLocalAssembler final
auto const flux = std::get<1>(flag_flux_dFlux);
auto const& dFlux = std::get<2>(flag_flux_dFlux);

auto const& wp = Base::_integration_method.getWeightedPoint(ip);
auto const w = sm.detJ * wp.getWeight() * sm.integralMeasure;
local_rhs.noalias() += sm.N * (flux * w);
local_rhs.noalias() += N * (flux * w);

if (static_cast<int>(dFlux.size()) != num_comp_total)
{
@@ -170,7 +167,7 @@ class PythonBoundaryConditionLocalAssembler final
// The assignement -= takes into account the sign convention
// of 1st-order in time ODE systems in OpenGeoSys.
local_Jac.block(top, left, width, height).noalias() -=
sm.N.transpose() * (dFlux[comp] * w) * sm.N;
N.transpose() * (dFlux[comp] * w) * N;
}
}
}
@@ -193,7 +190,6 @@ class PythonBoundaryConditionLocalAssembler final

private:
PythonBoundaryConditionData const& _data;
MeshLib::Element const& _element;
};

} // namespace ProcessLib
@@ -54,22 +54,22 @@ class RobinBoundaryConditionLocalAssembler final
unsigned const n_integration_points =
Base::_integration_method.getNumberOfPoints();

SpatialPosition pos;
pos.setElementID(id);
typename Base::NodalVectorType const alpha =
_data.alpha.getNodalValuesOnElement(Base::_element, t);

This comment has been minimized.

Copy link
@TomFischer

TomFischer Feb 25, 2019

Member

Can alpha be const? ✔️

typename Base::NodalVectorType const u_0 =
_data.u_0.getNodalValuesOnElement(Base::_element, t);

for (unsigned ip = 0; ip < n_integration_points; ++ip)
{
pos.setIntegrationPoint(ip);
auto const& ip_data = Base::_ns_and_weights[ip];

double const alpha = _data.alpha(t, pos)[0];
double const u_0 = _data.u_0(t, pos)[0];
auto const& N = ip_data.N;
auto const& w = ip_data.weight;

// flux = alpha * ( u_0 - u )
// adding a alpha term to the diagonal of the stiffness matrix
// and a alpha * u_0 term to the rhs vector
_local_K.diagonal().noalias() += ip_data.N * alpha * ip_data.weight;
_local_rhs.noalias() += ip_data.N * alpha * u_0 * ip_data.weight;
_local_K.diagonal().noalias() += N * alpha.dot(N) * w;
_local_rhs.noalias() += N * alpha.dot(N) * u_0.dot(N) * w;
}

auto const indices = NumLib::getIndices(id, dof_table_boundary);
@@ -74,7 +74,8 @@ class VariableDependentNeumannBoundaryConditionLocalAssembler final
NodalVectorType const coefficient_mixed_variables_node_values =
_data.coefficient_mixed_variables.getNodalValuesOnElement(
Base::_element, t);
unsigned const n_integration_points = _local_matrix_size;
unsigned const n_integration_points =
Base::_integration_method.getNumberOfPoints();

auto const indices_current_variable =
NumLib::getIndices(mesh_item_id, dof_table_boundary);
ProTip! Use n and p to navigate between commits in a pull request.
You can’t perform that action at this time.