Skip to content

Commit

Permalink
reduction to one shape matrix N
Browse files Browse the repository at this point in the history 
clang-format reformatting of branch
  • Loading branch information
bathmann committed Sep 7, 2018
1 parent 9a16ea1 commit a6b8f86
Showing 1 changed file with 45 additions and 69 deletions.
114 changes: 45 additions & 69 deletions ProcessLib/ComponentTransport/ComponentTransportFEM.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,19 +31,14 @@ namespace ComponentTransport
template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
struct IntegrationPointData final
{
IntegrationPointData(NodalRowVectorType const& N_C_,
GlobalDimNodalMatrixType const& dNdx_C_,
NodalRowVectorType const& N_p_,
GlobalDimNodalMatrixType const& dNdx_p_,
double const& integration_weight_)
: N_C(N_C_), dNdx_C(dNdx_C_), N_p(N_p_), dNdx_p(dNdx_p_),
integration_weight(integration_weight_)
IntegrationPointData(NodalRowVectorType const& N_,
GlobalDimNodalMatrixType const& dNdx_,
double const& integration_weight_)
: N(N_), dNdx(dNdx_), integration_weight(integration_weight_)
{}

NodalRowVectorType const N_C;
GlobalDimNodalMatrixType const dNdx_C;
NodalRowVectorType const N_p;
GlobalDimNodalMatrixType const dNdx_p;
NodalRowVectorType const N;
GlobalDimNodalMatrixType const dNdx;
double const integration_weight;

EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
Expand Down Expand Up @@ -109,21 +104,10 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
IntegrationMethod, GlobalDim>(
element, is_axially_symmetric, _integration_method);

auto const shape_matrices_C =
initShapeMatrices<ShapeFunction, ShapeMatricesType,
IntegrationMethod, GlobalDim>(
element, is_axially_symmetric, _integration_method);

auto const shape_matrices_p =
initShapeMatrices<ShapeFunction, ShapeMatricesType,
IntegrationMethod, GlobalDim>(
element, is_axially_symmetric, _integration_method);

for (unsigned ip = 0; ip < n_integration_points; ip++)
{
_ip_data.emplace_back(
shape_matrices_C[ip].N, shape_matrices_C[ip].dNdx,
shape_matrices_p[ip].N, shape_matrices_p[ip].dNdx,
shape_matrices[ip].N, shape_matrices[ip].dNdx,
_integration_method.getWeightedPoint(ip).getWeight() *
shape_matrices[ip].integralMeasure *
shape_matrices[ip].detJ);
Expand Down Expand Up @@ -172,8 +156,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
local_K.template block<num_nodes, num_nodes>(num_nodes, num_nodes);
auto Mpp =
local_M.template block<num_nodes, num_nodes>(num_nodes, num_nodes);
auto MpC =
local_M.template block<num_nodes, num_nodes>(num_nodes, 0);
auto MpC = local_M.template block<num_nodes, num_nodes>(num_nodes, 0);
auto Bp = local_b.template block<num_nodes, 1>(num_nodes, 0);

for (std::size_t ip(0); ip < n_integration_points; ++ip)
Expand All @@ -182,18 +165,16 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface


auto const& ip_data = _ip_data[ip];
auto const& N_C = ip_data.N_C;
auto const& dNdx_C = ip_data.dNdx_C;
auto const& N_p = ip_data.N_p;
auto const& dNdx_p = ip_data.dNdx_p;
auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
auto const& w = ip_data.integration_weight;

double C_int_pt = 0.0;
double p_int_pt = 0.0;
// Order matters: First C, then p!

NumLib::shapeFunctionInterpolate(C_nodal_values, N_C, C_int_pt);
NumLib::shapeFunctionInterpolate(p_nodal_values, N_p, p_int_pt);
NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt, p_int_pt);

// porosity model
auto const porosity =
_process_data.porous_media_properties.getPorosity(t, pos)
Expand Down Expand Up @@ -229,17 +210,16 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
GlobalDimVectorType const velocity =
_process_data.has_gravity
? GlobalDimVectorType(-K_over_mu *
(dNdx_p * p_nodal_values - density * b))
: GlobalDimVectorType(-K_over_mu * dNdx_p * p_nodal_values);
(dNdx * p_nodal_values - density * b))
: GlobalDimVectorType(-K_over_mu * dNdx * p_nodal_values);


const double drho_dp = _process_data.fluid_properties->getdValue(
MaterialLib::Fluid::FluidPropertyType::Density,
MaterialLib::Fluid::FluidPropertyType::Density,
vars,
MaterialLib::Fluid::PropertyVariableType::p);

const double drho_dC = _process_data.fluid_properties->getdValue(
MaterialLib::Fluid::FluidPropertyType::Density,
MaterialLib::Fluid::FluidPropertyType::Density,
vars,
MaterialLib::Fluid::PropertyVariableType::C);
double const velocity_magnitude = velocity.norm();
Expand All @@ -248,7 +228,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
? GlobalDimMatrixType(
(porosity * molecular_diffusion_coefficient +
solute_dispersivity_transverse *
velocity_magnitude) *
velocity_magnitude) *
I +
(solute_dispersivity_longitudinal -
solute_dispersivity_transverse) /
Expand All @@ -261,28 +241,25 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
I);

// matrix assembly
MCp.noalias() +=
w * N_C.transpose() * N_C * C_nodal_values *
retardation_factor * porosity * drho_dp * N_p;
MCC.noalias() +=
w * N_C.transpose() * retardation_factor * porosity *
(density + C_int_pt * drho_dC) * N_C;
MCp.noalias() += w * N.transpose() * N * C_nodal_values *
retardation_factor * porosity * drho_dp * N;
MCC.noalias() += w * N.transpose() * retardation_factor * porosity *
density * N +
w * N.transpose() * retardation_factor * porosity *
N * C_nodal_values * drho_dC * N;

KCC.noalias() +=
(-dNdx_C.transpose() * velocity * density * N_C +
dNdx_C.transpose() * density * hydrodynamic_dispersion *
dNdx_C + N_C.transpose() * decay_rate * porosity *
retardation_factor *
density * N_C) *
(-dNdx.transpose() * velocity * density * N +
dNdx.transpose() * density * hydrodynamic_dispersion * dNdx +
N.transpose() * decay_rate * porosity * retardation_factor *
density * N) *
w;

Mpp.noalias() += w * N_p.transpose() * porosity * drho_dp * N_p;
MpC.noalias() += w * N_p.transpose() * porosity * drho_dC * N_C;
Kpp.noalias() += w * dNdx_p.transpose() * density * K_over_mu *
dNdx_p;

Mpp.noalias() += w * N.transpose() * porosity * drho_dp * N;
MpC.noalias() += w * N.transpose() * porosity * drho_dC * N;
Kpp.noalias() += w * dNdx.transpose() * density * K_over_mu * dNdx;
if (_process_data.has_gravity)
Bp += w * density * density * dNdx_p.transpose() * K_over_mu *
b;
Bp += w * density * density * dNdx.transpose() * K_over_mu * b;
/* with Oberbeck-Boussing assumption density difference only exists
* in buoyancy effects */
}
Expand Down Expand Up @@ -311,16 +288,14 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface

MaterialLib::Fluid::FluidProperty::ArrayType vars;
auto const num_nodes = ShapeFunction::NPOINTS;
auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
&local_x[num_nodes], num_nodes);
auto C_nodal_values =
Eigen::Map<const NodalVectorType>(&local_x[0], num_nodes);
auto const p_nodal_values =
Eigen::Map<const NodalVectorType>(&local_x[num_nodes], num_nodes);

for (unsigned ip = 0; ip < n_integration_points; ++ip)
{
auto const& ip_data = _ip_data[ip];
auto const& N_p = ip_data.N_p;
auto const& dNdx_p = ip_data.dNdx_p;
auto const& N_C = ip_data.N_C;
auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;

pos.setIntegrationPoint(ip);

Expand All @@ -331,14 +306,15 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
GlobalDimMatrixType const K_over_mu = K / mu;

cache_mat.col(ip).noalias() = -K_over_mu * dNdx_p * p_nodal_values;
cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
if (_process_data.has_gravity)
{
double C_int_pt = 0.0;
double p_int_pt = 0.0;

NumLib::shapeFunctionInterpolate(C_nodal_values, N_C, C_int_pt);
NumLib::shapeFunctionInterpolate(p_nodal_values, N_p, p_int_pt);
NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt,
p_int_pt);

vars[static_cast<int>(
MaterialLib::Fluid::PropertyVariableType::C)] = C_int_pt;
vars[static_cast<int>(
Expand All @@ -359,10 +335,10 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
auto const& N_C = _ip_data[integration_point].N_C;
auto const& N = _ip_data[integration_point].N;

// assumes N is stored contiguously in memory
return Eigen::Map<const Eigen::RowVectorXd>(N_C.data(), N_C.size());
return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
}

private:
Expand Down

0 comments on commit a6b8f86

Please sign in to comment.