NumLib::detail Namespace Reference

Detailed Description

Classes
struct	ByComponent
struct	ByGlobalIndex
struct	ByLocation
struct	ByLocationAndComponent
struct	ComputeDerivativeWrtOneScalar_CD
struct	ComputeDerivativeWrtOneScalar_FD
struct	DefaultPerturbationStrategy
struct	FieldType
struct	IsScalar
struct	IsScalar< Eigen::Matrix< double, N, 1, Eigen::ColMajor, N, 1 > >
struct	Line
struct	LineByLocationAndComponentComparator
struct	LineByLocationComparator
struct	LowerOrderShapeFunctionOrSame
struct	LowerOrderShapeFunctionOrSame< ShapeFunction, std::void_t< typename NumLib::LowerDim< ShapeFunction >::type > >
struct	NumberOfDofs
struct	NumberOfDofs< Vectorial< ShapeFunction, Dim > >
struct	ShapeDataFieldType
struct	ShapeFunctionTraits

Typedefs
using	ComponentGlobalIndexDict
template<typename ElementTraitsLagrange>
using	GetMeshElement = typename ElementTraitsLagrange::Element
template<typename ElementTraitsLagrange>
using	GetShapeFunctionHigherOrder
template<typename ElementTraitsLagrange>
using	GetShapeFunctionLowerOrder
template<typename ShapeFunction>
using	GetShapeMatrixPolicy
template<typename ShapeMatrixPolicy>
using	GetShapeMatrix_N = typename ShapeMatrixPolicy::ShapeMatrices::ShapeType

Functions
template<typename... Ns_t, typename ElementDOFVector, int... Offsets, int... Sizes>
auto	localDOFImpl (ElementDOFVector const &x, boost::mp11::mp_list_c< int, Offsets... >, boost::mp11::mp_list_c< int, Sizes... >)
template<class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
void	computeMappingMatrices (const MeshLib::Element &, const double *natural_pt, const MeshLib::ElementCoordinatesMappingLocal &, T_SHAPE_MATRICES &shapemat, FieldType< ShapeMatrixType::N >)
template<class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
void	computeMappingMatrices (const MeshLib::Element &, const double *natural_pt, const MeshLib::ElementCoordinatesMappingLocal &, T_SHAPE_MATRICES &shapemat, FieldType< ShapeMatrixType::DNDR >)
static void	checkJacobianDeterminant (const double detJ, MeshLib::Element const &element)
template<class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
void	computeMappingMatrices (const MeshLib::Element &ele, const double *natural_pt, const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord, T_SHAPE_MATRICES &shapemat, FieldType< ShapeMatrixType::DNDR_J >)
template<class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
void	computeMappingMatrices (const MeshLib::Element &ele, const double *natural_pt, const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord, T_SHAPE_MATRICES &shapemat, FieldType< ShapeMatrixType::N_J >)
template<class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
void	computeMappingMatrices (const MeshLib::Element &ele, const double *natural_pt, const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord, T_SHAPE_MATRICES &shapemat, FieldType< ShapeMatrixType::DNDX >)
template<class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
void	computeMappingMatrices (const MeshLib::Element &ele, const double *natural_pt, const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord, T_SHAPE_MATRICES &shapemat, FieldType< ShapeMatrixType::ALL >)
template<class T_SHAPE_FUNC, class T_SHAPE_MATRICES, ShapeMatrixType T_SHAPE_MATRIX_TYPE>
void	naturalCoordinatesMappingComputeShapeMatrices (const MeshLib::Element &ele, const double *natural_pt, T_SHAPE_MATRICES &shapemat, const unsigned global_dim)
	Used to explicitly instantiate the NaturalCoordinatesMapping class template.
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapeHex20)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapeHex8)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapeLine2)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapeLine3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapePoint1)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapePrism15)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapePrism6)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapePyra13)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapePyra5)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapeQuad4)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapeQuad8)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapeQuad9)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapeTet10)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapeTet4)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapeTri3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN (ShapeTri6)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapePoint1, 0)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeLine2, 1)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeLine3, 1)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapePoint1, 1)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeLine2, 2)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeLine3, 2)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapePoint1, 2)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeQuad4, 2)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeQuad8, 2)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeQuad9, 2)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeTri3, 2)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeTri6, 2)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeHex20, 3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeHex8, 3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeLine2, 3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeLine3, 3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapePoint1, 3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapePrism15, 3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapePrism6, 3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapePyra13, 3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapePyra5, 3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeQuad4, 3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeQuad8, 3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeQuad9, 3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeTet10, 3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeTet4, 3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeTri3, 3)
	OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX (ShapeTri6, 3)
template<class T>
void	setMatrixZero (T &mat)
template<class T>
void	setVectorZero (T &vec)
template<class T_N, class T_DNDR, class T_J, class T_DNDX>
void	setZero (ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX > &shape, ShapeDataFieldType< ShapeMatrixType::N >)
template<class T_N, class T_DNDR, class T_J, class T_DNDX>
void	setZero (ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX > &shape, ShapeDataFieldType< ShapeMatrixType::DNDR >)
template<class T_N, class T_DNDR, class T_J, class T_DNDX>
void	setZero (ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX > &shape, ShapeDataFieldType< ShapeMatrixType::DNDR_J >)
template<class T_N, class T_DNDR, class T_J, class T_DNDX>
void	setZero (ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX > &shape, ShapeDataFieldType< ShapeMatrixType::N_J >)
template<class T_N, class T_DNDR, class T_J, class T_DNDX>
void	setZero (ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX > &shape, ShapeDataFieldType< ShapeMatrixType::DNDX >)
template<class T_N, class T_DNDR, class T_J, class T_DNDX>
void	setZero (ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX > &shape, ShapeDataFieldType< ShapeMatrixType::ALL >)
template<unsigned DOFOffset, typename NodalValues, typename ShapeMatrix>
void	shapeFunctionInterpolate (const NodalValues &, const ShapeMatrix &)
template<unsigned DOFOffset, typename NodalValues, typename ShapeMatrix, typename... ScalarTypes>
void	shapeFunctionInterpolate (const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
template<std::size_t IndexInTuple, typename Tuple>
double	getScalarOrVectorComponent (Tuple const &tuple, Eigen::Index component)
template<typename MeshElementType, typename IPData, typename FluxVectorType, typename Derived>
void	assembleAdvectionMatrix (IPData const &ip_data_vector, NumLib::ShapeMatrixCache const &shape_matrix_cache, std::vector< FluxVectorType > const &ip_flux_vector, Eigen::MatrixBase< Derived > &laplacian_matrix)
template<typename IPData, typename FluxVectorType, typename Derived>
void	assembleAdvectionMatrix (IPData const &ip_data_vector, std::vector< FluxVectorType > const &ip_flux_vector, Eigen::MatrixBase< Derived > &laplacian_matrix)
template<typename Derived>
void	applyFullUpwind (Eigen::VectorXd const &quasi_nodal_flux, Eigen::MatrixBase< Derived > &laplacian_matrix)
template<typename IPData, typename FluxVectorType, typename Derived>
void	applyFullUpwind (IPData const &ip_data_vector, std::vector< FluxVectorType > const &ip_flux_vector, Eigen::MatrixBase< Derived > &laplacian_matrix)
template<typename MatrixVectorType>
std::unique_ptr< MatrixVectorType >	newZeroedInstance (MathLib::MatrixSpecifications const &matrix_specification)
void	finalize (Mat &M)
void	calculateFluxCorrectedTransportPETSc (const double t, const double dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, const MathLib::MatrixSpecifications &matrix_specification, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b)
void	calculateFluxCorrectedTransport (const double t, const double dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, const MathLib::MatrixSpecifications &matrix_specification, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b)
Eigen::MatrixXd	getHydrodynamicDispersion (Eigen::MatrixXd const &pore_diffusion_coefficient, Eigen::VectorXd const &velocity, double const porosity, double const solute_dispersivity_transverse, double const solute_dispersivity_longitudinal)
Eigen::MatrixXd	getHydrodynamicDispersionWithArtificialDiffusion (IsotropicDiffusionStabilization const &stabilizer, std::size_t const element_id, Eigen::MatrixXd const &pore_diffusion_coefficient, Eigen::VectorXd const &velocity, double const porosity, double const solute_dispersivity_transverse, double const solute_dispersivity_longitudinal)
Eigen::VectorXd	computeAndersonWeights (Eigen::MatrixXd G)
bool	solvePicard (GlobalLinearSolver &linear_solver, GlobalMatrix &A, GlobalVector &rhs, GlobalVector &x, MathLib::LinearSolverBehaviour const linear_solver_behaviour)

Typedef Documentation

ComponentGlobalIndexDict

using NumLib::detail::ComponentGlobalIndexDict

Initial value:

 boost::multi_index::multi_index_container<
    Line, boost::multi_index::indexed_by<
              boost::multi_index::ordered_unique<
                  boost::multi_index::tag<ByLocationAndComponent>,
                  boost::multi_index::identity<Line>,
                  LineByLocationAndComponentComparator>,
              boost::multi_index::ordered_non_unique<
                  boost::multi_index::tag<ByLocation>,
                  boost::multi_index::identity<Line>, LineByLocationComparator>,
              boost::multi_index::ordered_non_unique<
                  boost::multi_index::tag<ByComponent>,
                  boost::multi_index::member<Line, int, &Line::comp_id>>,
              boost::multi_index::ordered_non_unique<
                  boost::multi_index::tag<ByGlobalIndex>,
                  boost::multi_index::member<Line, GlobalIndexType,
                                             &Line::global_index>>>>

Definition at line 95 of file ComponentGlobalIndexDict.h.

GetMeshElement

template<typename ElementTraitsLagrange>

using NumLib::detail::GetMeshElement = typename ElementTraitsLagrange::Element

Definition at line 16 of file ShapeMatrixCache.h.

GetShapeFunctionHigherOrder

template<typename ElementTraitsLagrange>

using NumLib::detail::GetShapeFunctionHigherOrder

Initial value:

 
    typename ElementTraitsLagrange::ShapeFunction

Definition at line 19 of file ShapeMatrixCache.h.

GetShapeFunctionLowerOrder

template<typename ElementTraitsLagrange>

using NumLib::detail::GetShapeFunctionLowerOrder

Initial value:

 
    typename ElementTraitsLagrange::LowerOrderShapeFunction

Definition at line 23 of file ShapeMatrixCache.h.

GetShapeMatrix_N

template<typename ShapeMatrixPolicy>

using NumLib::detail::GetShapeMatrix_N = typename ShapeMatrixPolicy::ShapeMatrices::ShapeType

Definition at line 34 of file ShapeMatrixCache.h.

GetShapeMatrixPolicy

template<typename ShapeFunction>

using NumLib::detail::GetShapeMatrixPolicy

Initial value:

 
    
    
    
    ShapeMatrixPolicyType<ShapeFunction, 3 >

Definition at line 27 of file ShapeMatrixCache.h.

Function Documentation

applyFullUpwind() [1/2]

template<typename Derived>

void NumLib::detail::applyFullUpwind	(	Eigen::VectorXd const &	quasi_nodal_flux,
		Eigen::MatrixBase< Derived > &	laplacian_matrix )

Definition at line 62 of file AdvectionMatrixAssembler.h.

{
    Eigen::VectorXd const down_mask =
        (quasi_nodal_flux.array() < 0).cast<double>();
    Eigen::VectorXd const down = quasi_nodal_flux.cwiseProduct(down_mask);
 
    double const q_in = -down.sum();
    if (q_in < std::numeric_limits<double>::epsilon())
    {
        return;
    }
 
    Eigen::VectorXd const up_mask =
        (quasi_nodal_flux.array() >= 0).cast<double>();
    Eigen::VectorXd const up = quasi_nodal_flux.cwiseProduct(up_mask);
 
    laplacian_matrix.diagonal().noalias() += up;
    laplacian_matrix.noalias() += down * up.transpose() / q_in;
}

Referenced by applyFullUpwind(), NumLib::assembleAdvectionMatrix(), and NumLib::assembleAdvectionMatrix().

applyFullUpwind() [2/2]

template<typename IPData, typename FluxVectorType, typename Derived>

void NumLib::detail::applyFullUpwind	(	IPData const &	ip_data_vector,
		std::vector< FluxVectorType > const &	ip_flux_vector,
		Eigen::MatrixBase< Derived > &	laplacian_matrix )

Definition at line 84 of file AdvectionMatrixAssembler.h.

{
    Eigen::VectorXd quasi_nodal_flux =
        Eigen::VectorXd::Zero(laplacian_matrix.rows());
 
    for (std::size_t ip = 0; ip < ip_flux_vector.size(); ++ip)
    {
        auto const& ip_data = ip_data_vector[ip];
        auto const w = ip_data.integration_weight;
        auto const& dNdx = ip_data.dNdx;
 
        quasi_nodal_flux.noalias() -= ip_flux_vector[ip].transpose() * dNdx * w;
    }
 
    applyFullUpwind(quasi_nodal_flux, laplacian_matrix);
}

References applyFullUpwind().

assembleAdvectionMatrix() [1/2]

template<typename MeshElementType, typename IPData, typename FluxVectorType, typename Derived>

void NumLib::detail::assembleAdvectionMatrix	(	IPData const &	ip_data_vector,
		NumLib::ShapeMatrixCache const &	shape_matrix_cache,
		std::vector< FluxVectorType > const &	ip_flux_vector,
		Eigen::MatrixBase< Derived > &	laplacian_matrix )

Definition at line 22 of file AdvectionMatrixAssembler.h.

{
    auto const& Ns = shape_matrix_cache.NsHigherOrder<MeshElementType>();
 
    for (std::size_t ip = 0; ip < ip_flux_vector.size(); ++ip)
    {
        auto const& ip_data = ip_data_vector[ip];
        auto const w = ip_data.integration_weight;
        auto const& dNdx = ip_data.dNdx;
        auto const& N = Ns[ip];
        laplacian_matrix.noalias() +=
            (N.transpose() * ip_flux_vector[ip].transpose() * dNdx * w)
#if defined(__GNUC__) && (__GNUC__ == 14) && (__GNUC_MINOR__ == 2) && \
    (__GNUC_PATCHLEVEL__ == 1)
                .eval()
#endif
            ;
    }
}

References NumLib::N, and NumLib::ShapeMatrixCache::NsHigherOrder().

Referenced by ProcessLib::HT::MonolithicHTFEM< ShapeFunction, GlobalDim >::assemble(), NumLib::assembleAdvectionMatrix(), NumLib::assembleAdvectionMatrix(), ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleBlockMatrices(), ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleComponentTransportEquation(), ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleHeatTransportEquation(), ProcessLib::HT::StaggeredHTFEM< ShapeFunction, GlobalDim >::assembleHeatTransportEquation(), ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobian(), and ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleWithJacobianComponentTransportEquation().

assembleAdvectionMatrix() [2/2]

template<typename IPData, typename FluxVectorType, typename Derived>

void NumLib::detail::assembleAdvectionMatrix	(	IPData const &	ip_data_vector,
		std::vector< FluxVectorType > const &	ip_flux_vector,
		Eigen::MatrixBase< Derived > &	laplacian_matrix )

Definition at line 46 of file AdvectionMatrixAssembler.h.

{
    for (std::size_t ip = 0; ip < ip_flux_vector.size(); ++ip)
    {
        auto const& ip_data = ip_data_vector[ip];
        auto const w = ip_data.integration_weight;
        auto const& N = ip_data.N;
        auto const& dNdx = ip_data.dNdx;
        laplacian_matrix.noalias() +=
            N.transpose() * ip_flux_vector[ip].transpose() * dNdx * w;
    }
}

References NumLib::N.

calculateFluxCorrectedTransport()

void NumLib::detail::calculateFluxCorrectedTransport	(	const double	t,
		const double	dt,
		std::vector< GlobalVector * > const &	x,
		std::vector< GlobalVector * > const &	x_prev,
		int const	process_id,
		const MathLib::MatrixSpecifications &	matrix_specification,
		GlobalMatrix &	M,
		GlobalMatrix &	K,
		GlobalVector &	b )

Definition at line 523 of file FluxCorrectedTransport.cpp.

{
#ifndef USE_PETSC
    calculateFluxCorrectedTransportSerial(t, dt, x, x_prev, process_id,
                                          matrix_specification, M, K, b);
#endif  // end of ifndef USE_PETSC
 
#ifdef USE_PETSC
    calculateFluxCorrectedTransportPETSc(t, dt, x, x_prev, process_id,
                                         matrix_specification, M, K, b);
#endif  // end of ifdef USE_PETSC
}

References calculateFluxCorrectedTransportPETSc().

Referenced by NumLib::computeFluxCorrectedTransport().

calculateFluxCorrectedTransportPETSc()

void NumLib::detail::calculateFluxCorrectedTransportPETSc	(	const double	t,
		const double	dt,
		std::vector< GlobalVector * > const &	x,
		std::vector< GlobalVector * > const &	x_prev,
		int const	process_id,
		const MathLib::MatrixSpecifications &	matrix_specification,
		GlobalMatrix &	M,
		GlobalMatrix &	K,
		GlobalVector &	b )

Definition at line 202 of file FluxCorrectedTransport.cpp.

{
    BaseLib::MPI::Mpi mpi(PETSC_COMM_WORLD);
    if (mpi.size > 1)
    {
        OGS_FATAL(
            "The current parallel implementation of flux corrected transport "
            "supports only one MPI rank, but the simulation uses {} ranks",
            mpi.size);
    }
    Mat K_raw = K.getRawMatrix();
    Mat D;
    PetscCallAbort(PETSC_COMM_WORLD,
                   MatDuplicate(K_raw, MAT_DO_NOT_COPY_VALUES, &D));
    PetscCallAbort(PETSC_COMM_WORLD, MatZeroEntries(D));
    finalize(D);
    PetscInt start, end;
    PetscCallAbort(PETSC_COMM_WORLD,
                   MatGetOwnershipRange(K.getRawMatrix(), &start, &end));
 
    PetscCallAbort(PETSC_COMM_WORLD, MatScale(K_raw, -1.0));
    K.finalizeAssembly();
    // compute artificial diffusion operator D
    for (PetscInt i = start; i < end; ++i)
    {
        PetscInt number_of_columns;
        PetscInt const* columns;
        PetscScalar const* values;
 
        PetscCallAbort(PETSC_COMM_WORLD, MatGetRow(K_raw, i, &number_of_columns,
                                                   &columns, &values));
 
        auto const columns_span =
            std::span<PetscInt const>(columns, number_of_columns);
        auto const values_span =
            std::span<PetscScalar const>(values, number_of_columns);
        for (auto const [j, kij] : std::views::zip(columns_span, values_span))
        {
            if (i == j)  // skip diagonal entries
            {
                continue;
            }
            PetscScalar kji;
            // MatGetValue is not supported for all PETSc matrix types
            PetscCallAbort(PETSC_COMM_WORLD, MatGetValue(K_raw, j, i, &kji));
            PetscScalar dij = std::max({-kij, 0.0, -kji});
            PetscCallAbort(PETSC_COMM_WORLD,
                           MatSetValue(D, i, j, dij, INSERT_VALUES));
        }
        // MatRestoreRow() is required after MatGetRow()
        PetscCallAbort(
            PETSC_COMM_WORLD,
            MatRestoreRow(K_raw, i, &number_of_columns, &columns, &values));
    }
    finalize(D);
    auto row_sums = newZeroedInstance<GlobalVector>(matrix_specification);
    Vec row_sums_raw = row_sums->getRawVector();
    PetscCallAbort(PETSC_COMM_WORLD, MatGetRowSum(D, row_sums_raw));
    PetscCallAbort(PETSC_COMM_WORLD, VecScale(row_sums_raw, -1.0));
    PetscCallAbort(PETSC_COMM_WORLD, VecAssemblyBegin(row_sums_raw));
    PetscCallAbort(PETSC_COMM_WORLD, VecAssemblyEnd(row_sums_raw));
    PetscCallAbort(PETSC_COMM_WORLD,
                   MatDiagonalSet(D, row_sums_raw, ADD_VALUES));
    finalize(D);
    // need to access x and x_prev in the following loops
    x[process_id]->setLocalAccessibleVector();
    x_prev[process_id]->setLocalAccessibleVector();
 
    // compute F
    Mat F;
    PetscCallAbort(PETSC_COMM_WORLD,
                   MatDuplicate(K.getRawMatrix(), MAT_DO_NOT_COPY_VALUES, &F));
    PetscCallAbort(PETSC_COMM_WORLD, MatZeroEntries(F));
    for (PetscInt i = start; i < end; ++i)
    {
        PetscInt cols;
        const PetscInt* col_indices;
        const PetscScalar* values;
 
        PetscCallAbort(PETSC_COMM_WORLD,
                       MatGetRow(D, i, &cols, &col_indices, &values));
 
        auto const columns = std::span<PetscInt const>(col_indices, cols);
        auto const values_D = std::span<PetscScalar const>(values, cols);
 
        for (auto const [j, dij] : std::views::zip(columns, values_D))
        {
            PetscScalar const xi = x[process_id]->get(i);
            PetscScalar const xj = x[process_id]->get(j);
            PetscScalar const fij = -dij * (xj - xi);
            PetscCallAbort(PETSC_COMM_WORLD,
                           MatSetValue(F, i, j, fij, INSERT_VALUES));
        }
        PetscCallAbort(PETSC_COMM_WORLD,
                       MatRestoreRow(D, i, &cols, &col_indices, &values));
    }
    finalize(F);
    auto& M_raw = M.getRawMatrix();
 
    auto xdot = newZeroedInstance<GlobalVector>(matrix_specification);
    xdot->setLocalAccessibleVector();
    for (PetscInt i = start; i < end; ++i)
    {
        xdot->add(i, (x[process_id]->get(i) - x_prev[process_id]->get(i)) / dt);
    }
    xdot->finalizeAssembly();
    xdot->setLocalAccessibleVector();
 
    for (PetscInt i = start; i < end; ++i)
    {
        PetscInt cols;
        PetscInt const* col_indices;
        PetscScalar const* values;
 
        PetscCallAbort(PETSC_COMM_WORLD,
                       MatGetRow(M_raw, i, &cols, &col_indices, &values));
 
        auto const columns = std::span<PetscInt const>(col_indices, cols);
        auto const values_M = std::span<PetscScalar const>(values, cols);
 
        for (auto const [j, mij] : std::views::zip(columns, values_M))
        {
            PetscScalar const fij = -mij * (xdot->get(j) - xdot->get(i));
            PetscCallAbort(PETSC_COMM_WORLD,
                           MatSetValue(F, i, j, fij, ADD_VALUES));
        }
        PetscCallAbort(PETSC_COMM_WORLD,
                       MatRestoreRow(M_raw, i, &cols, &col_indices, &values));
    }
    finalize(F);
    auto P_plus = newZeroedInstance<GlobalVector>(matrix_specification);
    auto P_minus = newZeroedInstance<GlobalVector>(matrix_specification);
 
    P_plus->setLocalAccessibleVector();
    P_minus->setLocalAccessibleVector();
    // Fill P_plus, P_minus, Q_plus, Q_minus
    for (PetscInt i = start; i < end; ++i)
    {
        PetscInt cols;
        PetscInt const* col_indices;
        PetscScalar const* values;
 
        PetscCallAbort(PETSC_COMM_WORLD,
                       MatGetRow(F, i, &cols, &col_indices, &values));
 
        auto const column_indices =
            std::span<PetscInt const>(col_indices, cols);
        auto const values_span = std::span<PetscScalar const>(values, cols);
        for (auto const [j, fij] : std::views::zip(column_indices, values_span))
        {
            // skip diagonal
            if (i == j)
            {
                continue;
            }
            P_plus->add(i, std::max(0., fij));
            P_minus->add(i, std::min(0., fij));
        }
        PetscCallAbort(PETSC_COMM_WORLD,
                       MatRestoreRow(F, i, &cols, &col_indices, &values));
    }
    P_plus->finalizeAssembly();
    P_plus->setLocalAccessibleVector();
    P_minus->finalizeAssembly();
    P_minus->setLocalAccessibleVector();
    auto Q_plus = newZeroedInstance<GlobalVector>(matrix_specification);
    auto Q_minus = newZeroedInstance<GlobalVector>(matrix_specification);
    Q_plus->setLocalAccessibleVector();
    Q_minus->setLocalAccessibleVector();
 
    for (PetscInt i = start; i < end; ++i)
    {
        PetscInt cols;
        PetscInt const* col_indices;
 
        PetscScalar Q_plus_i = 0.0;
        PetscScalar Q_minus_i = 0.0;
 
        PetscCallAbort(PETSC_COMM_WORLD,
                       MatGetRow(F, i, &cols, &col_indices, nullptr));
        auto const column_indices =
            std::span<PetscInt const>(col_indices, cols);
        for (auto const j : column_indices)
        {
            // skip diagonal
            if (i == j)
            {
                continue;
            }
            PetscScalar const x_prev_i = x_prev[process_id]->get(i);
            PetscScalar const x_prev_j = x_prev[process_id]->get(j);
 
            Q_plus_i = std::max({Q_plus_i, 0., x_prev_j - x_prev_i});
            Q_minus_i = std::min({Q_minus_i, 0., x_prev_j - x_prev_i});
        }
        Q_plus->set(i, Q_plus_i);
        Q_minus->set(i, Q_minus_i);
        PetscCallAbort(PETSC_COMM_WORLD,
                       MatRestoreRow(F, i, &cols, &col_indices, nullptr));
    }
    Q_plus->finalizeAssembly();
    Q_plus->setLocalAccessibleVector();
    Q_minus->finalizeAssembly();
    Q_minus->setLocalAccessibleVector();
 
    auto R_plus = newZeroedInstance<GlobalVector>(matrix_specification);
    R_plus->setLocalAccessibleVector();
 
    auto row_sums_M = newZeroedInstance<GlobalVector>(matrix_specification);
    Vec M_L = row_sums_M->getRawVector();
    PetscCallAbort(PETSC_COMM_WORLD, MatGetRowSum(M_raw, M_L));
    row_sums_M->finalizeAssembly();
    row_sums_M->setLocalAccessibleVector();
 
    for (auto k = R_plus->getRangeBegin(); k < R_plus->getRangeEnd(); ++k)
    {
        PetscScalar const P_plus_i = P_plus->get(k);
        if (P_plus_i == 0.0)
        {
            continue;
        }
        PetscScalar const mi = row_sums_M->get(k);
        PetscScalar const Q_plus_i = Q_plus->get(k);
        R_plus->set(k, std::min(1.0, mi * Q_plus_i / dt / P_plus_i));
    }
    R_plus->finalizeAssembly();
    R_plus->setLocalAccessibleVector();
    auto R_minus = newZeroedInstance<GlobalVector>(matrix_specification);
    R_minus->setLocalAccessibleVector();
 
    for (auto k = R_minus->getRangeBegin(); k < R_minus->getRangeEnd(); ++k)
    {
        PetscScalar const P_minus_i = P_minus->get(k);
        if (P_minus_i == 0.0)
        {
            continue;
        }
        PetscScalar const mi = row_sums_M->get(k);
        PetscScalar const Q_minus_i = Q_minus->get(k);
        R_minus->set(k, std::min(1.0, mi * Q_minus_i / dt / P_minus_i));
    }
    R_minus->finalizeAssembly();
    R_minus->setLocalAccessibleVector();
    // walk over F, R_plus, and R_minus and compute alpha values; set entries in
    // the rhs b that limit the antidiffusive flux
    for (PetscInt i = start; i < end; ++i)
    {
        PetscInt cols;
        PetscInt const* col_indices;
        PetscScalar const* values;
 
        PetscCallAbort(PETSC_COMM_WORLD,
                       MatGetRow(F, i, &cols, &col_indices, &values));
 
        auto col_indices_span = std::span<PetscInt const>(col_indices, cols);
        auto values_span = std::span<PetscScalar const>(values, cols);
 
        double alpha_ij;
 
        for (auto [j, fij] : std::views::zip(col_indices_span, values_span))
        {
            if (i == j)
            {
                continue;
            }
            if (fij > 0.)
            {
                alpha_ij = std::min(R_plus->get(i), R_minus->get(j));
            }
            else
            {
                alpha_ij = std::min(R_minus->get(i), R_plus->get(j));
            }
            b.add(i, alpha_ij * fij);
        }
        PetscCallAbort(PETSC_COMM_WORLD,
                       MatRestoreRow(F, i, &cols, &col_indices, &values));
    }
 
    // compute low-order operator
    // update K_raw with D: K_raw += D;
    for (PetscInt i = start; i < end; ++i)
    {
        PetscInt cols_D;
        PetscInt const* col_indices_D;
        PetscScalar const* values_D;
        PetscCallAbort(PETSC_COMM_WORLD,
                       MatGetRow(D, i, &cols_D, &col_indices_D, &values_D));
 
        auto const column_indices =
            std::span<PetscInt const>(col_indices_D, cols_D);
        auto const values = std::span<PetscScalar const>(values_D, cols_D);
        for (auto [j, Dij] : std::views::zip(column_indices, values))
        {
            PetscCallAbort(PETSC_COMM_WORLD,
                           MatSetValues(K_raw, 1, &i, 1, &j, &Dij, ADD_VALUES));
        }
        PetscCallAbort(PETSC_COMM_WORLD,
                       MatRestoreRow(D, i, &cols_D, &col_indices_D, &values_D));
    }
    finalize(K_raw);
    PetscCallAbort(PETSC_COMM_WORLD, MatScale(K_raw, -1.0));
    K.finalizeAssembly();
    M.setZero();
    for (PetscInt i = start; i < end; ++i)
    {
        PetscScalar const row_sum_M_i = row_sums_M->get(i);
        MatSetValues(M.getRawMatrix(), 1, &i, 1, &i, &row_sum_M_i,
                     INSERT_VALUES);
    }
    M.finalizeAssembly();
    MatDestroy(&D);
    MatDestroy(&F);
}

References MathLib::EigenVector::add(), finalize(), MathLib::EigenMatrix::getRawMatrix(), newZeroedInstance(), OGS_FATAL, MathLib::EigenMatrix::setZero(), and BaseLib::MPI::Mpi::size.

Referenced by calculateFluxCorrectedTransport().

checkJacobianDeterminant()

void NumLib::detail::checkJacobianDeterminant ( const double detJ, MeshLib::Element const & element )

static

Definition at line 85 of file NaturalCoordinatesMapping.cpp.

{
    if (detJ > 0)
    {  // The usual case
        return;
    }
 
    if (detJ < 0)
    {
        ERR("det J = {:g} is negative for element {:d}.",
            detJ,
            element.getID());
#ifndef NDEBUG
        std::cerr << element << "\n";
#endif  // NDEBUG
        OGS_FATAL(
            "Please check whether the node numbering of the element is correct,"
            "or additional elements (like boundary elements) are still present "
            "in the mesh.");
    }
 
    if (detJ == 0)
    {
        ERR("det J is zero for element {:d}.", element.getID());
#ifndef NDEBUG
        std::cerr << element << "\n";
#endif  // NDEBUG
        OGS_FATAL(
            "Please check whether:\n"
            "\t the element nodes may have the same coordinates,\n"
            "\t or the coordinates of all nodes are not given on the x-axis "
            "for a 1D problem or in the xy-plane in the 2D case.\n"
            "The first case can occur, if not all boundary elements"
            "were removed from the bulk mesh.");
    }
}

References ERR(), MeshLib::Element::getID(), and OGS_FATAL.

Referenced by computeMappingMatrices(), and computeMappingMatrices().

computeAndersonWeights()

Eigen::VectorXd NumLib::detail::computeAndersonWeights

(

Eigen::MatrixXd

G

)

Computes the Anderson mixing weights \( \theta \) minimising \( \|\sum_i \theta_i f_i\|^2 \) subject to \( \sum_i \theta_i = 1 \), given the Gram matrix \( G = F^T F \) of the stored residuals \( f_i = g(x_i) - x_i \).

The constrained least-squares problem is solved via the Lagrange formulation

\[ \begin{bmatrix} G & \mathbf{1} \\ \mathbf{1}^T & 0 \end{bmatrix} \begin{bmatrix} \theta \\ \lambda \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ 1 \end{bmatrix} \]

G is normalized internally by its largest diagonal entry ( \( \max_i \|f_i\|^2 \)) to keep the augmented system well-conditioned regardless of the residuals' magnitude; the solution \( \theta \) is invariant under this scaling, because the scale factor cancels in \( \theta = G^{-1}\mathbf{1} / (\mathbf{1}^T G^{-1}\mathbf{1}) \).

Linearly dependent (or nearly dependent) residuals are the classical failure mode of Anderson acceleration. The computed mixture is therefore accepted only if

1. it predicts a smaller residual norm than the plain step it would replace, which rules out the arbitrary minimizer that a degenerate history admits, and
2. its weights stay bounded, which rules out the mixtures of near-identical iterates whose value is pure cancellation error, and
3. its predicted residual is not numerically zero while every stored step still has a non-negligible norm, which rules out the collinear steps of differing magnitude whose exact cancellation needs only modest weights and therefore escapes the two tests above.

Otherwise \( \theta = (0,\dots,0,1) \) is returned, i.e. unit weight on the newest stored step, which reproduces the plain (damped) Picard update. The same fallback applies when all stored residuals vanish. Acceleration thus degrades to plain Picard instead of amplifying rounding error.

Note

This is an implementation detail of AndersonAcceleration, declared in a separate header so that it stays out of the public class header while remaining reachable from the unit tests.

Parameters

G	taken by value on purpose: it is normalized in place (divided by its largest diagonal entry) before the solve, so the caller's matrix must not be touched.

Precondition

G is symmetric positive semi-definite and at least 1x1.

Definition at line 21 of file AndersonAcceleration.cpp.

{
    int const history_size = static_cast<int>(G.rows());
 
    // At least one residual is required: G is the Gram matrix of the stored
    // steps, so maxCoeff() below is a reduction over a non-empty diagonal.
    if (history_size < 1)
    {
        OGS_FATAL(
            "Anderson acceleration: the mixing weights were requested for an "
            "empty history (Gram matrix of size {:d}x{:d}).",
            history_size, static_cast<int>(G.cols()));
    }
 
    // Unit weight on the newest stored step, i.e. the plain (damped) Picard
    // update. Used whenever the stored history does not admit trustworthy
    // mixing weights.
    auto const newest_step_only = [history_size]()
    {
        Eigen::VectorXd theta = Eigen::VectorXd::Zero(history_size);
        theta(history_size - 1) = 1.0;
        return theta;
    };
 
    double const g_scale = G.diagonal().maxCoeff();
    // Negated comparison, so that a NaN scale takes this branch as well.
    if (!(g_scale > 0.0))
    {
        // All stored steps vanish; there is nothing left to mix.
        return newest_step_only();
    }
    G /= g_scale;
 
    Eigen::MatrixXd M(history_size + 1, history_size + 1);
    M.topLeftCorner(history_size, history_size) = G;
    M.topRightCorner(history_size, 1).setOnes();
    M.bottomLeftCorner(1, history_size).setOnes();
    M(history_size, history_size) = 0.0;
 
    Eigen::VectorXd rhs_aa(history_size + 1);
    rhs_aa.setZero();
    rhs_aa(history_size) = 1.0;
 
    Eigen::VectorXd const theta =
        M.fullPivLu().solve(rhs_aa).head(history_size);
 
    if (!theta.allFinite())
    {
        // A fallback silently changes the iterate the solver would otherwise
        // take, so it is reported at INFO level rather than hidden in DBUG.
        INFO(
            "Anderson acceleration: the mixing weights came out non-finite. "
            "Falling back to the plain Picard step for this iteration.");
        return newest_step_only();
    }
 
    // The rank-revealing solve above never fails outright, so the weights have
    // to be validated on their own merits. Note that an invertibility test on
    // M would be the wrong check: it rejects healthy histories whose residual
    // norms span many orders of magnitude (the normal situation for a
    // converging iteration) while passing the ill-conditioned cases that
    // actually do harm.
    //
    // (a) Descent: the mixture only earns its place if the residual norm it
    //     predicts is smaller than that of the plain step it would replace.
    //     For a degenerate history - linearly dependent steps, in particular
    //     duplicates - the minimizer is not unique and the solve returns an
    //     arbitrary one, which this test discards.
    //     Both norms are read off the same normalized G (divided by g_scale
    //     above), so the shared 1/g_scale factor cancels and the comparison is
    //     exactly the one on the unscaled residuals.
    double const mixed_residual_norm_2 = theta.dot(G * theta);
    double const newest_step_norm_2 = G(history_size - 1, history_size - 1);
 
    // (b) No long lever arms: weights far outside [0, 1] mean the mixed
    //     iterate is a difference of near-identical vectors, i.e. dominated by
    //     cancellation. Healthy histories stay at |theta| ~ 1, whereas steps
    //     that agree to k digits produce weights of magnitude 10^k.
    //
    //     The threshold is a heuristic, not a derived bound. A weight of
    //     magnitude 10^k sacrifices about k of the ~16 significant decimal
    //     digits of a double to cancellation; capping at 10^2 admits the
    //     modest lever arms of genuinely useful mixing (empirically |theta| up
    //     to ~10) while rejecting the 10^3-and-up weights that signal a
    //     degenerate history. It is deliberately loose: the descent test (a)
    //     is the primary guard, and this one only catches the cancellation
    //     cases that slip past it.
    constexpr double max_weight = 1e2;
 
    // (c) No exact cancellation of substantial steps: guards (a) and (b) both
    //     miss the case of exactly (or near-exactly) collinear stored steps
    //     whose magnitudes differ enough that cancelling them needs only
    //     modest weights - e.g. two steps in a 2x ratio need theta = (2, -1),
    //     well inside the max_weight cap, yet a rank-deficient G lets their
    //     combination be driven to an exactly-zero model residual regardless
    //     of how far either step actually is from the fixed point.
    //     That is trustworthy only when it is explained by a stored step that
    //     is already that small on its own (the ZeroResidualGetsFullWeight
    //     case, where G has a genuine zero diagonal entry): if every stored
    //     step still has a non-negligible norm, an exactly-zero mixed
    //     residual can only be an algebraic artefact of the steps happening
    //     to be parallel, not evidence of proximity to the fixed point.
    constexpr double numerically_zero = 1e-8;
    bool const exact_cancellation_of_substantial_steps =
        mixed_residual_norm_2 < numerically_zero &&
        G.diagonal().minCoeff() > numerically_zero;
 
    // Negated comparison, so that a NaN norm takes the fallback as well.
    if (!(mixed_residual_norm_2 < newest_step_norm_2) ||
        theta.cwiseAbs().maxCoeff() > max_weight ||
        exact_cancellation_of_substantial_steps)
    {
        // A fallback silently changes the iterate the solver would otherwise
        // take, so it is reported at INFO level rather than hidden in DBUG.
        INFO(
            "Anderson acceleration: rejected the mixture of {:d} stored steps "
            "(predicted residual {:g} vs. {:g} for the plain step, largest "
            "weight {:g}). Falling back to the plain Picard step for this "
            "iteration.",
            history_size, mixed_residual_norm_2, newest_step_norm_2,
            theta.cwiseAbs().maxCoeff());
        return newest_step_only();
    }
 
    return theta;
}

References INFO(), and OGS_FATAL.

Referenced by NumLib::AndersonAcceleration::accelerate().

computeMappingMatrices() [1/6]

template<class T_SHAPE_FUNC, class T_SHAPE_MATRICES>

void NumLib::detail::computeMappingMatrices ( const MeshLib::Element & , const double * natural_pt, const MeshLib::ElementCoordinatesMappingLocal & , T_SHAPE_MATRICES & shapemat, FieldType< ShapeMatrixType::DNDR >  )

inline

Definition at line 71 of file NaturalCoordinatesMapping.cpp.

{
    if constexpr (T_SHAPE_FUNC::DIM != 0)
    {
        double* const dNdr = shapemat.dNdr.data();
        T_SHAPE_FUNC::computeGradShapeFunction(natural_pt, dNdr);
    }
}

computeMappingMatrices() [2/6]

template<class T_SHAPE_FUNC, class T_SHAPE_MATRICES>

void NumLib::detail::computeMappingMatrices ( const MeshLib::Element & , const double * natural_pt, const MeshLib::ElementCoordinatesMappingLocal & , T_SHAPE_MATRICES & shapemat, FieldType< ShapeMatrixType::N >  )

inline

Definition at line 60 of file NaturalCoordinatesMapping.cpp.

{
    T_SHAPE_FUNC::computeShapeFunction(natural_pt, shapemat.N);
}

Referenced by computeMappingMatrices(), computeMappingMatrices(), computeMappingMatrices(), computeMappingMatrices(), and naturalCoordinatesMappingComputeShapeMatrices().

computeMappingMatrices() [3/6]

template<class T_SHAPE_FUNC, class T_SHAPE_MATRICES>

void NumLib::detail::computeMappingMatrices ( const MeshLib::Element & ele, const double * natural_pt, const MeshLib::ElementCoordinatesMappingLocal & ele_local_coord, T_SHAPE_MATRICES & shapemat, FieldType< ShapeMatrixType::ALL >  )

inline

Definition at line 229 of file NaturalCoordinatesMapping.cpp.

{
    computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
        ele,
        natural_pt,
        ele_local_coord,
        shapemat,
        FieldType<ShapeMatrixType::N>());
    computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
        ele,
        natural_pt,
        ele_local_coord,
        shapemat,
        FieldType<ShapeMatrixType::DNDX>());
}

References computeMappingMatrices().

computeMappingMatrices() [4/6]

template<class T_SHAPE_FUNC, class T_SHAPE_MATRICES>

void NumLib::detail::computeMappingMatrices ( const MeshLib::Element & ele, const double * natural_pt, const MeshLib::ElementCoordinatesMappingLocal & ele_local_coord, T_SHAPE_MATRICES & shapemat, FieldType< ShapeMatrixType::DNDR_J >  )

inline

Definition at line 124 of file NaturalCoordinatesMapping.cpp.

{
    if constexpr (T_SHAPE_FUNC::DIM != 0)
    {
        computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
            ele,
            natural_pt,
            ele_local_coord,
            shapemat,
            FieldType<ShapeMatrixType::DNDR>());
 
        auto constexpr nnodes = T_SHAPE_FUNC::NPOINTS;
 
        // jacobian: J=[dx/dr dy/dr // dx/ds dy/ds]
        for (auto k = decltype(nnodes){0}; k < nnodes; k++)
        {
            const MathLib::Point3d& mapped_pt =
                ele_local_coord.getMappedCoordinates(k);
            shapemat.J +=
                shapemat.dNdr.col(k) *
                mapped_pt.asEigenVector3d().head(T_SHAPE_FUNC::DIM).transpose();
        }
 
        shapemat.detJ = shapemat.J.determinant();
        checkJacobianDeterminant(shapemat.detJ, ele);
    }
    else
    {
        shapemat.detJ = 1.0;
    }
}

References MathLib::Point3d::asEigenVector3d(), checkJacobianDeterminant(), computeMappingMatrices(), and MeshLib::ElementCoordinatesMappingLocal::getMappedCoordinates().

computeMappingMatrices() [5/6]

template<class T_SHAPE_FUNC, class T_SHAPE_MATRICES>

void NumLib::detail::computeMappingMatrices ( const MeshLib::Element & ele, const double * natural_pt, const MeshLib::ElementCoordinatesMappingLocal & ele_local_coord, T_SHAPE_MATRICES & shapemat, FieldType< ShapeMatrixType::DNDX >  )

inline

Definition at line 184 of file NaturalCoordinatesMapping.cpp.

{
    computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
        ele,
        natural_pt,
        ele_local_coord,
        shapemat,
        FieldType<ShapeMatrixType::DNDR_J>());
    if constexpr (T_SHAPE_FUNC::DIM != 0)
    {
        checkJacobianDeterminant(shapemat.detJ, ele);
 
        // J^-1, dshape/dx
        shapemat.invJ.noalias() = shapemat.J.inverse();
 
        assert(shapemat.dNdr.rows() == ele.getDimension());
        const unsigned global_dim = ele_local_coord.getGlobalDimension();
        if (global_dim == T_SHAPE_FUNC::DIM)
        {
            shapemat.dNdx
                .template topLeftCorner<T_SHAPE_FUNC::DIM,
                                        T_SHAPE_FUNC::NPOINTS>()
                .noalias() = shapemat.invJ * shapemat.dNdr;
        }
        else
        {
            auto const& matR =
                (ele_local_coord.getRotationMatrixToGlobal().topLeftCorner(
                     global_dim, T_SHAPE_FUNC::DIM))
                    .eval();
 
            auto const invJ_dNdr = shapemat.invJ * shapemat.dNdr;
            auto const dshape_global = matR * invJ_dNdr;
            shapemat.dNdx =
                dshape_global.topLeftCorner(global_dim, T_SHAPE_FUNC::NPOINTS);
        }
    }
}

References checkJacobianDeterminant(), computeMappingMatrices(), MeshLib::Element::getDimension(), MeshLib::ElementCoordinatesMappingLocal::getGlobalDimension(), and MeshLib::ElementCoordinatesMappingLocal::getRotationMatrixToGlobal().

computeMappingMatrices() [6/6]

template<class T_SHAPE_FUNC, class T_SHAPE_MATRICES>

void NumLib::detail::computeMappingMatrices ( const MeshLib::Element & ele, const double * natural_pt, const MeshLib::ElementCoordinatesMappingLocal & ele_local_coord, T_SHAPE_MATRICES & shapemat, FieldType< ShapeMatrixType::N_J >  )

inline

Definition at line 162 of file NaturalCoordinatesMapping.cpp.

{
    computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
        ele,
        natural_pt,
        ele_local_coord,
        shapemat,
        FieldType<ShapeMatrixType::N>());
    computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
        ele,
        natural_pt,
        ele_local_coord,
        shapemat,
        FieldType<ShapeMatrixType::DNDR_J>());
}

References computeMappingMatrices().

finalize()

void NumLib::detail::finalize

(

Mat &

M

)

Definition at line 194 of file FluxCorrectedTransport.cpp.

{
    PetscCallAbort(PETSC_COMM_WORLD,
                   MatAssemblyBegin(M, MatAssemblyType::MAT_FINAL_ASSEMBLY));
    PetscCallAbort(PETSC_COMM_WORLD,
                   MatAssemblyEnd(M, MatAssemblyType::MAT_FINAL_ASSEMBLY));
}

Referenced by calculateFluxCorrectedTransportPETSc().

getHydrodynamicDispersion()

Eigen::MatrixXd NumLib::detail::getHydrodynamicDispersion ( Eigen::MatrixXd const & pore_diffusion_coefficient, Eigen::VectorXd const & velocity, double const porosity, double const solute_dispersivity_transverse, double const solute_dispersivity_longitudinal )

inline

Definition at line 15 of file HydrodynamicDispersion.h.

{
    double const velocity_magnitude = velocity.norm();
    if (velocity_magnitude == 0.0)
    {
        return porosity * pore_diffusion_coefficient;
    }
 
    auto const dim = velocity.size();
    Eigen::MatrixXd const& I(Eigen::MatrixXd::Identity(dim, dim));
    return porosity * pore_diffusion_coefficient +
           solute_dispersivity_transverse * velocity_magnitude * I +
           (solute_dispersivity_longitudinal - solute_dispersivity_transverse) /
               velocity_magnitude * velocity * velocity.transpose();
}

Referenced by NumLib::computeHydrodynamicDispersion().

getHydrodynamicDispersionWithArtificialDiffusion()

Eigen::MatrixXd NumLib::detail::getHydrodynamicDispersionWithArtificialDiffusion ( IsotropicDiffusionStabilization const & stabilizer, std::size_t const element_id, Eigen::MatrixXd const & pore_diffusion_coefficient, Eigen::VectorXd const & velocity, double const porosity, double const solute_dispersivity_transverse, double const solute_dispersivity_longitudinal )

inline

Definition at line 36 of file HydrodynamicDispersion.h.

{
    double const velocity_magnitude = velocity.norm();
    if (velocity_magnitude == 0.0)
    {
        return porosity * pore_diffusion_coefficient;
    }
 
    double const artificial_diffusion =
        stabilizer.computeArtificialDiffusion(element_id, velocity_magnitude);
 
    auto const dim = velocity.size();
    Eigen::MatrixXd const& I(Eigen::MatrixXd::Identity(dim, dim));
    return porosity * pore_diffusion_coefficient +
           (solute_dispersivity_transverse * velocity_magnitude +
            artificial_diffusion) *
               I +
           (solute_dispersivity_longitudinal - solute_dispersivity_transverse) /
               velocity_magnitude * velocity * velocity.transpose();
}

References NumLib::IsotropicDiffusionStabilization::computeArtificialDiffusion().

Referenced by NumLib::computeHydrodynamicDispersion().

getScalarOrVectorComponent()

template<std::size_t IndexInTuple, typename Tuple>

double NumLib::detail::getScalarOrVectorComponent	(	Tuple const &	tuple,
		Eigen::Index	component )

Definition at line 32 of file NumericalDifferentiation.h.

{
    auto const& value = std::get<IndexInTuple>(tuple);
 
    if constexpr (IsScalar<std::remove_cvref_t<decltype(value)>>::value)
    {
        return value;
    }
    else
    {
        return value[component];
    }
}

Referenced by NumLib::detail::ComputeDerivativeWrtOneScalar_CD::operator()(), and NumLib::detail::ComputeDerivativeWrtOneScalar_FD< Value >::operator()().

localDOFImpl()

template<typename... Ns_t, typename ElementDOFVector, int... Offsets, int... Sizes>

auto NumLib::detail::localDOFImpl	(	ElementDOFVector const &	x,
		boost::mp11::mp_list_c< int, Offsets... >	,
		boost::mp11::mp_list_c< int, Sizes... >	)

Definition at line 34 of file LocalDOF.h.

{
    static_assert(((Sizes > 0) && ...));
    static_assert(((Offsets >= 0) && ...));
 
    return std::tuple<Eigen::Map<typename MatrixPolicyType::VectorType<
        NumberOfDofs<Ns_t>::value> const>...>{{x.data() + Offsets, Sizes}...};
}

Referenced by NumLib::localDOF().

naturalCoordinatesMappingComputeShapeMatrices()

template<class T_SHAPE_FUNC, class T_SHAPE_MATRICES, ShapeMatrixType T_SHAPE_MATRIX_TYPE>

void NumLib::detail::naturalCoordinatesMappingComputeShapeMatrices	(	const MeshLib::Element &	ele,
		const double *	natural_pt,
		T_SHAPE_MATRICES &	shapemat,
		const unsigned	global_dim )

Used to explicitly instantiate the NaturalCoordinatesMapping class template.

Definition at line 253 of file NaturalCoordinatesMapping.cpp.

{
    const MeshLib::ElementCoordinatesMappingLocal ele_local_coord(ele,
                                                                  global_dim);
 
    detail::computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
        ele,
        natural_pt,
        ele_local_coord,
        shapemat,
        detail::FieldType<T_SHAPE_MATRIX_TYPE>());
}

References computeMappingMatrices().

Referenced by NumLib::NaturalCoordinatesMapping< ShapeFunctionType, ShapeMatrices >::computeShapeMatrices().

newZeroedInstance()

template<typename MatrixVectorType>

std::unique_ptr< MatrixVectorType > NumLib::detail::newZeroedInstance

(

MathLib::MatrixSpecifications const &

matrix_specification

)

Definition at line 24 of file FluxCorrectedTransport.cpp.

{
    auto result = MathLib::MatrixVectorTraits<MatrixVectorType>::newInstance(
        matrix_specification);
    result->setZero();
    return result;
}

Referenced by calculateFluxCorrectedTransportPETSc().

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [1/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapeHex20

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [2/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapeHex8

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [3/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapeLine2

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [4/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapeLine3

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [5/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapePoint1

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [6/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapePrism15

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [7/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapePrism6

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [8/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapePyra13

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [9/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapePyra5

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [10/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapeQuad4

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [11/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapeQuad8

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [12/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapeQuad9

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [13/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapeTet10

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [14/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapeTet4

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [15/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapeTri3

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN() [16/16]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN

(

ShapeTri6

)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [1/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeHex20	,
		3	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [2/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeHex8	,
		3	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [3/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeLine2	,
		1	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [4/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeLine2	,
		2	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [5/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeLine2	,
		3	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [6/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeLine3	,
		1	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [7/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeLine3	,
		2	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [8/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeLine3	,
		3	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [9/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapePoint1	,
		0	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [10/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapePoint1	,
		1	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [11/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapePoint1	,
		2	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [12/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapePoint1	,
		3	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [13/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapePrism15	,
		3	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [14/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapePrism6	,
		3	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [15/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapePyra13	,
		3	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [16/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapePyra5	,
		3	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [17/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeQuad4	,
		2	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [18/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeQuad4	,
		3	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [19/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeQuad8	,
		2	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [20/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeQuad8	,
		3	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [21/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeQuad9	,
		2	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [22/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeQuad9	,
		3	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [23/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeTet10	,
		3	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [24/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeTet4	,
		3	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [25/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeTri3	,
		2	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [26/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeTri3	,
		3	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [27/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeTri6	,
		2	)

OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX() [28/28]

NumLib::detail::OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX	(	ShapeTri6	,
		3	)

setMatrixZero()

template<class T>

void NumLib::detail::setMatrixZero

(

T &

mat

)

Definition at line 12 of file ShapeMatrices-impl.h.

{
    // mat.setZero();
    const std::size_t n = mat.rows() * mat.cols();
    auto* v = mat.data();
    for (std::size_t i = 0; i < n; i++)
    {
        v[i] = .0;
    }
}

Referenced by setZero(), setZero(), and setZero().

setVectorZero()

template<class T>

void NumLib::detail::setVectorZero

(

T &

vec

)

Definition at line 24 of file ShapeMatrices-impl.h.

{
    // vec.setZero();
    const std::size_t n = vec.size();
    auto* v = vec.data();
    for (std::size_t i = 0; i < n; i++)
    {
        v[i] = .0;
    }
}

Referenced by setZero().

setZero() [1/6]

template<class T_N, class T_DNDR, class T_J, class T_DNDX>

void NumLib::detail::setZero ( ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX > & shape, ShapeDataFieldType< ShapeMatrixType::ALL >  )

inline

Definition at line 86 of file ShapeMatrices-impl.h.

{
    setZero(shape, ShapeDataFieldType<ShapeMatrixType::N>());
    setZero(shape, ShapeDataFieldType<ShapeMatrixType::DNDX>());
}

References setZero().

setZero() [2/6]

template<class T_N, class T_DNDR, class T_J, class T_DNDX>

void NumLib::detail::setZero ( ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX > & shape, ShapeDataFieldType< ShapeMatrixType::DNDR >  )

inline

Definition at line 52 of file ShapeMatrices-impl.h.

{
    setMatrixZero(shape.dNdr);
}

References NumLib::ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX >::dNdr, and setMatrixZero().

setZero() [3/6]

template<class T_N, class T_DNDR, class T_J, class T_DNDX>

void NumLib::detail::setZero ( ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX > & shape, ShapeDataFieldType< ShapeMatrixType::DNDR_J >  )

inline

Definition at line 59 of file ShapeMatrices-impl.h.

{
    setZero(shape, ShapeDataFieldType<ShapeMatrixType::DNDR>());
    setMatrixZero(shape.J);
    shape.detJ = .0;
    shape.integralMeasure = 0.0;
}

References NumLib::ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX >::detJ, NumLib::ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX >::integralMeasure, NumLib::ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX >::J, setMatrixZero(), and setZero().

setZero() [4/6]

template<class T_N, class T_DNDR, class T_J, class T_DNDX>

void NumLib::detail::setZero ( ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX > & shape, ShapeDataFieldType< ShapeMatrixType::DNDX >  )

inline

Definition at line 77 of file ShapeMatrices-impl.h.

{
    setZero(shape, ShapeDataFieldType<ShapeMatrixType::DNDR_J>());
    setMatrixZero(shape.invJ);
    setMatrixZero(shape.dNdx);
}

References NumLib::ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX >::dNdx, NumLib::ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX >::invJ, setMatrixZero(), and setZero().

setZero() [5/6]

template<class T_N, class T_DNDR, class T_J, class T_DNDX>

void NumLib::detail::setZero ( ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX > & shape, ShapeDataFieldType< ShapeMatrixType::N >  )

inline

Definition at line 45 of file ShapeMatrices-impl.h.

{
    setVectorZero(shape.N);
}

References NumLib::ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX >::N, and setVectorZero().

Referenced by setZero(), setZero(), setZero(), setZero(), and NumLib::ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX >::setZero().

setZero() [6/6]

template<class T_N, class T_DNDR, class T_J, class T_DNDX>

void NumLib::detail::setZero ( ShapeMatrices< T_N, T_DNDR, T_J, T_DNDX > & shape, ShapeDataFieldType< ShapeMatrixType::N_J >  )

inline

Definition at line 69 of file ShapeMatrices-impl.h.

{
    setZero(shape, ShapeDataFieldType<ShapeMatrixType::N>());
    setZero(shape, ShapeDataFieldType<ShapeMatrixType::DNDR_J>());
}

References setZero().

shapeFunctionInterpolate() [1/2]

template<unsigned DOFOffset, typename NodalValues, typename ShapeMatrix>

void NumLib::detail::shapeFunctionInterpolate	(	const NodalValues &	,
		const ShapeMatrix &	)

See also

NumLib::shapeFunctionInterpolate()

Definition at line 20 of file Interpolation.h.

{
}

Referenced by ProcessLib::HCNonAdvectiveFreeComponentFlowBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >::assemble(), ProcessLib::HeatConduction::LocalAssemblerData< ShapeFunction, GlobalDim >::assemble(), ProcessLib::HeatTransportBHE::HeatTransportBHELocalAssemblerSoil< ShapeFunction >::assemble(), ProcessLib::HT::MonolithicHTFEM< ShapeFunction, GlobalDim >::assemble(), ProcessLib::RichardsComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assemble(), ProcessLib::RichardsFlow::LocalAssemblerData< ShapeFunction, GlobalDim >::assemble(), ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assemble(), ProcessLib::SteadyStateDiffusion::LocalAssemblerData< ShapeFunction, GlobalDim >::assemble(), ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::assemble(), ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::assemble(), ProcessLib::TwoPhaseFlowWithPP::TwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::assemble(), ProcessLib::VariableDependentNeumannBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >::assemble(), ProcessLib::WellboreCompensateNeumannBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >::assemble(), ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::assemble(), ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleBlockMatrices(), ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleHeatTransportEquation(), ProcessLib::HT::StaggeredHTFEM< ShapeFunction, GlobalDim >::assembleHeatTransportEquation(), ProcessLib::HT::StaggeredHTFEM< ShapeFunction, GlobalDim >::assembleHydraulicEquation(), ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleReactionEquationConcrete(), ProcessLib::HeatConduction::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleWithJacobian(), ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobian(), ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::assembleWithJacobian(), ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForDeformationEquations(), ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::calculateIntPtDarcyVelocity(), ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::calculateIntPtLiquidDensity(), ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::computeSecondaryVariableConcrete(), ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::computeSecondaryVariableConcrete(), ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::computeSecondaryVariableConcrete(), ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::getFlux(), ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >::getFlux(), ProcessLib::LiquidFlow::LiquidFlowLocalAssembler< ShapeFunction, GlobalDim >::getFlux(), ProcessLib::SteadyStateDiffusion::LocalAssemblerData< ShapeFunction, GlobalDim >::getFlux(), ProcessLib::RichardsComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::getIntPtDarcyVelocity(), ProcessLib::RichardsFlow::LocalAssemblerData< ShapeFunction, GlobalDim >::getIntPtDarcyVelocity(), ProcessLib::SteadyStateDiffusion::LocalAssemblerData< ShapeFunction, GlobalDim >::getIntPtDarcyVelocity(), ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >::getIntPtDarcyVelocityLocal(), ProcessLib::RichardsComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::getIntPtSaturation(), ProcessLib::EmbeddedAnchor< GlobalDim >::getShapeMatricesAndGlobalIndicesAndDisplacements(), ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::initializeChemicalSystemConcrete(), ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::setChemicalSystemConcrete(), ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::setInitialConditionsConcrete(), ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::setInitialConditionsConcrete(), ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, DisplacementDim, ConstitutiveTraits >::setInitialConditionsConcrete(), shapeFunctionInterpolate(), and NumLib::shapeFunctionInterpolate().

shapeFunctionInterpolate() [2/2]

template<unsigned DOFOffset, typename NodalValues, typename ShapeMatrix, typename... ScalarTypes>

void NumLib::detail::shapeFunctionInterpolate	(	const NodalValues &	nodal_values,
		const ShapeMatrix &	shape_matrix_N,
		double &	interpolated_value,
		ScalarTypes &...	interpolated_values )

See also

NumLib::shapeFunctionInterpolate()

Definition at line 28 of file Interpolation.h.

{
    auto const num_nodes = shape_matrix_N.size();
    double iv = 0.0;
 
    for (auto n = decltype(num_nodes){0}; n < num_nodes; ++n)
    {
        iv += nodal_values[DOFOffset * num_nodes + n] * shape_matrix_N[n];
    }
 
    interpolated_value = iv;
 
    shapeFunctionInterpolate<DOFOffset + 1>(nodal_values, shape_matrix_N,
                                            interpolated_values...);
}

References shapeFunctionInterpolate().

solvePicard()

bool NumLib::detail::solvePicard	(	GlobalLinearSolver &	linear_solver,
		GlobalMatrix &	A,
		GlobalVector &	rhs,
		GlobalVector &	x,
		MathLib::LinearSolverBehaviour const	linear_solver_behaviour )

Definition at line 51 of file NonlinearSolver.cpp.

{
    if (linear_solver_behaviour ==
            MathLib::LinearSolverBehaviour::RECOMPUTE_AND_STORE ||
        linear_solver_behaviour == MathLib::LinearSolverBehaviour::REUSE)
    {
        WARN(
            "The performance optimization to skip the linear solver compute() "
            "step is not implemented for PETSc or LIS linear solvers.");
    }
 
    BaseLib::RunTime time_linear_solver;
    time_linear_solver.start();
 
    bool const iteration_succeeded = linear_solver.solve(A, rhs, x);
 
    INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
 
    if (iteration_succeeded)
    {
        return true;
    }
 
    ERR("Picard: The linear solver failed in the solve() step.");
    return false;
}

References BaseLib::RunTime::elapsed(), ERR(), INFO(), MathLib::RECOMPUTE_AND_STORE, MathLib::REUSE, MathLib::EigenLisLinearSolver::solve(), BaseLib::RunTime::start(), and WARN().

Referenced by NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::solve().