#include <ImplicitVectorOperators_fwd.hpp>

Detailed Description

template<class shape_storage_type, class matrix_type, class rhs_type>
class mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >

This class represents implicit vector operators that fill given matrix M and right hand side rhs with appropriate coefficients approximating differential operators with shape functions from given shape storage ss.

They are intended for implicit solutions of vector PDEs. Component specific equations can also be set using eq.

Warning: The operators grad(), value() etc. are lazy and will not write to the matrix until summed with another operator or assigned to. If eager evaluation is required, use the eval() method on each operation.; Each evaluated operation is only added to the matrix and rhs. Setting the same equation twice thus has the effect of multiplying it by 2. Therefore, the matrix and the rhs should be initialized (usually to zero) before setting any equations.

Template Parameters

shape_storage_type	Any shape storage class satisfying the Shape storage concept.
matrix_type	Usually and Eigen sparse matrix.
rhs_type	Usually an Eigen `VectorXd`.

The matrix and rhs should be of size at least N*dim, where N is the number of nodes in the domain and dim its dimensionality.

See also: ImplicitOperators

Usage example:

    Eigen::SparseMatrix<double, Eigen::RowMajor> M(2 * N, 2 * N);
    M.reserve(Range<int>(2 * N, 2 * support_size));
    Eigen::VectorXd field(2 * N), rhs(2 * N);
 
    auto op = storage.implicitVectorOperators(M, rhs);
    std::cout << op << std::endl;
 
    // Set equation on interior
    for (int i : domain.interior()) {
        double x = domain.pos(i, 0), y = domain.pos(i, 1);
        op.lap(i) + op.graddiv(i) = {2+4*x*y-8*PI*PI*std::sin(2*PI*x), 2*(2*x*x + y*y)};
    }
    // Set boundary conditions
    for (int i : domain.types() == -1) {
        double y = domain.pos(i, 1);
        op.value(i) = {y*y, 0};
    }
    for (int i : domain.types() == -2) {
        double y = domain.pos(i, 1);
        op.value(i) = {y*y, y*y};
    }
    for (int i : domain.types() == -3) {
        double x = domain.pos(i, 0);
        op.value(i) = {std::sin(2*PI*x), 0};
    }
    for (int i : domain.types() == -4) {
        double x = domain.pos(i, 0);
        op.value(i) = {1+std::sin(2*PI*x), x*x};
    }

Definition at line 44 of file ImplicitVectorOperators_fwd.hpp.

Public Member Functions
	ImplicitVectorOperators ()
	Default constructor sets offset to `0` and pointers to `nullptr`. More...

	ImplicitVectorOperators (const shape_storage_t &ss)
	Set only shape storage, the rest as default constructor. More...

	ImplicitVectorOperators (const shape_storage_t &ss, matrix_t &M, rhs_t &rhs, int row_offset=0, int col_offset=0)
	Set shape storage and problem matrix and rhs. More...

void	setProblem (matrix_t &M, rhs_t &rhs, int row_offset=0, int col_offset=0)
	Sets current matrix and right hand side. More...

void	setRowOffset (int row_offset)
	Sets row offset for given matrix, treating it as is the first row had index `row_offset`. More...

void	setColOffset (int col_offset)
	Sets col offset for given matrix, treating it as is the first column had index `col_offset`. More...

bool	hasShapes () const
	Returns `true` if operators have a non-null pointer to storage and `false` otherwise. More...

bool	hasMatrix () const
	Returns `true` if operators have a non-null pointer to problem matrix. More...

bool	hasRhs () const
	Returns `true` if operators have a non-null pointer to problem right hand side. More...

Equation	eq (int num)
	Choose one specific component of the vector equation to write to. More...

ValueOp	value (int node, int row)
	Sets implicit equation that value of a vector field is equal to some other value. More...

ValueOp	value (int node)
	Same as value with row equal to current node. More...

BasicOp< Lap< dim > >	lap (int node, int row)
	Sets implicit equation that Laplacian of a vector field is equal to some other value. More...

BasicOp< Lap< dim > >	lap (int node)
	Same as lap with row equal to current node. More...

template<typename op_family_t >
BasicOp< op_family_t >	apply (int node, typename op_family_t::operator_t o, int row)
	Add the weights for operator `o` to the appropriate elements in the matrix. More...

template<typename op_family_t >
BasicOp< op_family_t >	apply (int node, typename op_family_t::operator_t o)
	Same as apply with row equal to current node. More...

template<typename op_family_t >
BasicOp< op_family_t >	apply (int node)
	Overload for default-constructible operator. More...

GradOp	grad (int node, const domain_vector_t &v, int row)
	Sets implicit equation that gradient of a vector field along `v` is equal to some other value. More...

GradOp	grad (int node, const domain_vector_t &v)
	Same as grad with row equal to current node. More...

GradDivOp	graddiv (int node, int row)
	Sets implicit equation that gradient of divergence of a vector field is equal to some other value. More...

GradDivOp	graddiv (int node)
	Same as graddiv with row equal to current node. More...

GradOp	neumann (int node, const domain_vector_t &unit_normal, int row)
	Sets neumann boundary conditions in node `node`. More...

GradOp	neumann (int node, const domain_vector_t &n)
	Same as neumann with row equal to current node. More...

TractionOp	traction (int node, scalar_t lam, scalar_t mu, const domain_vector_t &unit_normal, int row)
	Sets traction boundary conditions in node `node`. More...

TractionOp	traction (int node, scalar_t lam, scalar_t mu, const domain_vector_t &n)
	Same as traction with row equal to current node. More...

Public Types
enum	{ dim = shape_storage_t::dim }
	Store dimension of the domain. More...

typedef shape_storage_type	shape_storage_t
	Type of shape storage. More...

typedef matrix_type	matrix_t
	Matrix type. More...

typedef rhs_type	rhs_t
	Right hand side type. More...

typedef shape_storage_t::vector_t	domain_vector_t
	Vector type for vectors in the domain, usually float or double. More...

typedef matrix_t::Scalar	scalar_t
	Scalar type of matrix elements. More...

typedef Eigen::Matrix< scalar_t, dim, 1 >	vector_t
	Vector type. More...

Classes
class	BasicOp
	Class representing basic operators. More...

class	Equation
	Represents one scalar component of a vector equation. More...

class	GradDivOp
	Class representing the gradient of the divergence operator, i.e. operator \(\nabla\nabla\cdot\). More...

class	GradOp
	Class representing the directional derivative (gradient) operation. More...

class	RowVecOp
	Class representing an operation on a specific set of rows of the matrix. More...

class	TractionOp
	Class representing the traction operator. Useful for setting boundary conditions. More...

class	ValueOp
	Class representing the "evaluate" operation, i.e. the zero-th derivative. More...

class	VecOpBase
	Base class for all elementary implicit vector operations. More...

Friends
template<typename S , typename M , typename R >
std::ostream &	operator<< (std::ostream &os, const ImplicitVectorOperators< S, M, R > &op)
	Output basic info about given operators. More...

Private Member Functions
void	addToM (int i, int j, scalar_t v)
	Adds `v` to `M(i, j)`, i.e. `M(i, j) += v`. More...

void	addToRhs (int i, const vector_t &v)
	Adds components of `v` to `i`, `i+N`, ..., `i+(dim-1)*N`-th rows of `rhs`. More...

Private Attributes
const shape_storage_t *	ss
	Shape storage, but name is shortened for readability. More...

matrix_t *	M
	Pointer to problem matrix. More...

rhs_t *	rhs
	Pointer to right hand side. More...

int	row_offset
	Row offset to be used when accessing matrix or rhs coefficients. More...

int	col_offset
	Column offset to be used when accessing matrix coefficients. More...

Member Enumeration Documentation

◆ anonymous enum

template<class shape_storage_type , class matrix_type , class rhs_type >

anonymous enum

Store dimension of the domain.

Enumerator
dim	Dimensionality of the domain.

Definition at line 55 of file ImplicitVectorOperators_fwd.hpp.

Constructor & Destructor Documentation

◆ ImplicitVectorOperators() [1/3]

template<class shape_storage_type , class matrix_type , class rhs_type >

mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::ImplicitVectorOperators ( )

inline

Default constructor sets offset to 0 and pointers to nullptr.

Definition at line 319 of file ImplicitVectorOperators_fwd.hpp.

◆ ImplicitVectorOperators() [2/3]

template<class shape_storage_type , class matrix_type , class rhs_type >

mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::ImplicitVectorOperators ( const shape_storage_t & ss )

inlineexplicit

Set only shape storage, the rest as default constructor.

Definition at line 322 of file ImplicitVectorOperators_fwd.hpp.

◆ ImplicitVectorOperators() [3/3]

template<class shape_storage_type , class matrix_type , class rhs_type >

mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::ImplicitVectorOperators	(	const shape_storage_t &	ss,
		matrix_t &	M,
		rhs_t &	rhs,
		int	row_offset = `0`,
		int	col_offset = `0`
	)

Set shape storage and problem matrix and rhs.

Parameters

ss	Class storing all computed shapes.
M	Problem matrix. Must have at least `ss->size()*dim` rows and cols.
rhs	Problem right hand side. Must have at least `ss->size()*dim` rows.
row_offset	Instead of counting rows from 0, count them from row `row_offset`.
col_offset	Instead of counting columns from 0, count them from row `col_offset`.

Warning: Matrices M and rhs have values only added to them and should be zero initialized by the user. Since this is a common mistake, a warning is printed if this is not the case when in debug mode.

This class is usually constructed directly from shape storage using the implicitVectorOperators() member function.

Definition at line 16 of file ImplicitVectorOperators.hpp.

Member Function Documentation

◆ addToM()

template<class shape_storage_type , class matrix_type , class rhs_type >

void mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::addToM	(	int	i,
		int	j,
		scalar_t	v
	)

inlineprivate

Adds v to M(i, j), i.e. M(i, j) += v.

Definition at line 68 of file ImplicitVectorOperators_fwd.hpp.

◆ addToRhs()

template<class shape_storage_type , class matrix_type , class rhs_type >

void mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::addToRhs	(	int	i,
		const vector_t &	v
	)

inlineprivate

Adds components of v to i, i+N, ..., i+(dim-1)*N-th rows of rhs.

Definition at line 74 of file ImplicitVectorOperators_fwd.hpp.

◆ apply() [1/3]

template<class shape_storage_type , class matrix_type , class rhs_type >

template<typename op_family_t >

BasicOp<op_family_t> mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::apply ( int node )

inline

Overload for default-constructible operator.

Definition at line 445 of file ImplicitVectorOperators_fwd.hpp.

◆ apply() [2/3]

template<class shape_storage_type , class matrix_type , class rhs_type >

template<typename op_family_t >

BasicOp<op_family_t> mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::apply	(	int	node,
		typename op_family_t::operator_t	o
	)

inline

Same as apply with row equal to current node.

Definition at line 439 of file ImplicitVectorOperators_fwd.hpp.

◆ apply() [3/3]

template<class shape_storage_type , class matrix_type , class rhs_type >

template<typename op_family_t >

BasicOp<op_family_t> mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::apply	(	int	node,
		typename op_family_t::operator_t	o,
		int	row
	)

inline

Add the weights for operator o to the appropriate elements in the matrix.

Parameters

node
o	Operator family
row	Write in this matrix row.

Returns: This function is lazy and returns an "operation" which will write the weights in the matrix, when evaluated.

Definition at line 433 of file ImplicitVectorOperators_fwd.hpp.

◆ eq()

template<class shape_storage_type , class matrix_type , class rhs_type >

Equation mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::eq ( int num )

inline

Choose one specific component of the vector equation to write to.

Usage example:

alpha*op.eq(0).c(0).der1(i, 0) + beta*op.eq(0).c(1).der1(i, 1) = tx(domain.pos(i));

sets equation \(\alpha \dpar{u}{x}(i) + \beta \dpar{v}{y}(i) = t_x(i)\) as the first scalar equation in the matrix.

See also: Equation::c

Definition at line 374 of file ImplicitVectorOperators_fwd.hpp.

◆ grad() [1/2]

template<class shape_storage_type , class matrix_type , class rhs_type >

GradOp mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::grad	(	int	node,
		const domain_vector_t &	v
	)

inline

Same as grad with row equal to current node.

Definition at line 467 of file ImplicitVectorOperators_fwd.hpp.

◆ grad() [2/2]

template<class shape_storage_type , class matrix_type , class rhs_type >

GradOp mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::grad	(	int	node,
		const domain_vector_t &	v,
		int	row
	)

inline

Sets implicit equation that gradient of a vector field along v is equal to some other value.

The code alpha*op.grad(node, v) = r fills the node-th, node+N-th and node+2*N-th row of matrix M these rows of equation \( M u = r \) are a good approximation of the equation

\[ \alpha (\vec{v} \cdot \nabla) \vec{u}(p) = \vec{r}(p), \]

where \(p\) is the node-th point and \(\vec{v}\) is the vector v.

Note: Above example was made for dim = 3, other dimensions are analogous.

Parameters

node	Index of a node from 0 to `N` for which to write the equation for. This means only `node`-th`, node+N`-th and `node+2*N`-th rows of the matrix are changed. User must make sure that the matrix is large enough.
row	Write equation in this specific row. Row with index `node` is chosen by default.
v	Vector to multiply the gradient with.

Example:

    // Solve problem for u unknown vector field:
    // grad u . v + 1/2 laplace u = f
    // grad u . (y, x) + 1/2 laplace u = (x^3+y+2 x y^2, x^2+(y-1) y)
    // on [0, 1] x [0, 1] with BC:
    // u(x, 0) = (0, -x)
    // u(x, 1) = (x*x, 0)
    // u(0, y) = (0, 0)
    // u(1, y) = (y, y-1)
    // with solution:
    // u(x, y) = (x^2 y, yx - x)
    std::function<Vec2d(Vec2d)> analytical = [] (const Vec2d& p) {
        return Vec2d(p[0]*p[0]*p[1], (p[1]-1)*p[0]);
    };
    // Prepare domain
    BoxShape<Vec2d> box(0., 1.);
    DomainDiscretization<Vec2d> domain = box.discretizeWithStep(0.05);
    int N = domain.size();
    int support_size = 9;
    domain.findSupport(FindClosest(support_size));
    // Prepare operators and matrix
    WLS<Monomials<Vec2d>, GaussianWeight<Vec2d>> wls(2, 15.0 / N);
    auto storage = domain.computeShapes(wls);
    Eigen::SparseMatrix<double> M(2*N, 2*N);
    M.reserve(storage.supportSizes()+storage.supportSizes());
    Eigen::VectorXd rhs(2*N); rhs.setZero();
    auto op = storage.implicitVectorOperators(M, rhs);
 
    Range<Vec2d> v(N);
    for (int i = 0; i < N; ++i) {
        v[i] = Vec2d({domain.pos(i, 1), domain.pos(i, 0)});
    }
    // Set equation on interior
    for (int i : domain.interior()) {
        double x = domain.pos(i, 0), y = domain.pos(i, 1);
        op.grad(i, v[i]) + 0.5*op.lap(i) = {x*x*x + y + 2*x*y*y, x*x + (y-1)*y};
    }
    // Set boundary conditions
    for (int i : domain.types() == -1) {
        op.value(i) = 0;
    }
    for (int i : domain.types() == -2) {
        double y = domain.pos(i, 1);
        op.value(i) = {y, y-1};
    }
    for (int i : domain.types() == -3) {
        double x = domain.pos(i, 0);
        op.value(i) = {0, -x};
    }
    for (int i : domain.types() == -4) {
        double x = domain.pos(i, 0);
        op.value(i) = {x*x, 0};
    }
 
    // Sparse solve
    Eigen::BiCGSTAB<Eigen::SparseMatrix<double, Eigen::RowMajor>> solver;
    solver.compute(M);
    Eigen::VectorXd sol = solver.solve(rhs);
    for (int i = 0; i < N; ++i) {
        Vec2d correct = analytical(domain.pos(i));
        ASSERT_NEAR(sol[i], correct[0], 1e-3);
        ASSERT_NEAR(sol[i+N], correct[1], 1e-3);
    }

Definition at line 465 of file ImplicitVectorOperators_fwd.hpp.

◆ graddiv() [1/2]

template<class shape_storage_type , class matrix_type , class rhs_type >

GradDivOp mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::graddiv ( int node )

inline

Same as graddiv with row equal to current node.

Definition at line 488 of file ImplicitVectorOperators_fwd.hpp.

◆ graddiv() [2/2]

template<class shape_storage_type , class matrix_type , class rhs_type >

GradDivOp mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::graddiv	(	int	node,
		int	row
	)

inline

Sets implicit equation that gradient of divergence of a vector field is equal to some other value.

The code alpha*op.graddiv(node) = v fills the node-th, node+N-th and node+2*N-th row of matrix M these rows of equation \( M u = v \) are a good approximation of the equation

\[ \alpha \nabla(\nabla \cdot \vec{u})(p) = \vec{v}(p), \]

where \(p\) is the node-th point.

Note: Above example was made for dim = 3, other dimensions are analogous.

Parameters

node	Index of a node from 0 to `N` for which to write the equation for. This means only `node`-th`, node+N`-th and `node+2*N`-th rows of the matrix are changed. User must make sure that the matrix is large enough.
row	Write equation in this specific row. Row with index `node` is chosen by default.

Example:

    // Solve problem for u unknown vector field:
    // grad div u + laplace u = (2 + 4xy - 8pi^2sin(2*pi*x), 2(2x^2+y^2))
    // on [0, 1] x [0, 1] with BC:
    // u(0, y) = (y^2, 0)
    // u(1, y) = (y^2, y^2)
    // u(x, 0) = (sin(2*pi*x), 0)
    // u(x, 1) = (1+sin(2*pi*x, x^2)
    // with solution:
    // u(x, y) = (y^2 + sin(2*pi*x), x^2y^2)
    std::function<Vec2d(Vec2d)> analytical = [] (const Vec2d& p) {
        return Vec2d({p[1]*p[1]+std::sin(2*PI*p[0]), p[0]*p[0]*p[1]*p[1]});
    };
    // Prepare domain
    BoxShape<Vec2d> box(0., 1.);
    DomainDiscretization<Vec2d> domain = box.discretizeWithStep(0.025);
    int N = domain.size();
    int support_size = 6;
    domain.findSupport(FindClosest(support_size));
    // Prepare operators and matrix
    WLS<Monomials<Vec2d>, GaussianWeight<Vec2d>> wls(2, std::sqrt(15.0 / N));
    auto storage = domain.computeShapes<sh::graddiv|sh::lap>(wls, domain.interior());
 
    Eigen::SparseMatrix<double, Eigen::RowMajor> M(2*N, 2*N);
    M.reserve(Range<int>(2*N, support_size));
    Eigen::VectorXd rhs(2*N); rhs.setZero();
 
    auto op = storage.implicitVectorOperators(M, rhs);
    // Set equation on interior
    for (int i : domain.interior()) {
        double x = domain.pos(i, 0), y = domain.pos(i, 1);
        op.lap(i) + op.graddiv(i) = {2+4*x*y-8*PI*PI*std::sin(2*PI*x), 2*(2*x*x + y*y)};
    }
    // Set boundary conditions
    for (int i : domain.types() == -1) {
        double y = domain.pos(i, 1);
        op.value(i) = {y*y, 0};
    }
    for (int i : domain.types() == -2) {
        double y = domain.pos(i, 1);
        op.value(i) = {y*y, y*y};
    }
    for (int i : domain.types() == -3) {
        double x = domain.pos(i, 0);
        op.value(i) = {std::sin(2*PI*x), 0};
    }
    for (int i : domain.types() == -4) {
        double x = domain.pos(i, 0);
        op.value(i) = {1+std::sin(2*PI*x), x*x};
    }
 
    // Sparse solve
    Eigen::BiCGSTAB<Eigen::SparseMatrix<double, Eigen::RowMajor>> solver;
    solver.compute(M);
    Eigen::VectorXd sol = solver.solve(rhs);
    for (int i = 0; i < N; ++i) {
        Vec2d correct = analytical(domain.pos(i));
        EXPECT_NEAR(sol[i], correct[0], 0.8e-2);
        EXPECT_NEAR(sol[i+N], correct[1], 0.8e-2);
    }

Definition at line 486 of file ImplicitVectorOperators_fwd.hpp.

◆ hasMatrix()

template<class shape_storage_type , class matrix_type , class rhs_type >

bool mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::hasMatrix ( ) const

inline

Returns true if operators have a non-null pointer to problem matrix.

Definition at line 363 of file ImplicitVectorOperators_fwd.hpp.

◆ hasRhs()

template<class shape_storage_type , class matrix_type , class rhs_type >

bool mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::hasRhs ( ) const

inline

Returns true if operators have a non-null pointer to problem right hand side.

Definition at line 365 of file ImplicitVectorOperators_fwd.hpp.

◆ hasShapes()

template<class shape_storage_type , class matrix_type , class rhs_type >

bool mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::hasShapes ( ) const

inline

Returns true if operators have a non-null pointer to storage and false otherwise.

Definition at line 361 of file ImplicitVectorOperators_fwd.hpp.

◆ lap() [1/2]

template<class shape_storage_type , class matrix_type , class rhs_type >

BasicOp<Lap<dim> > mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::lap ( int node )

inline

Same as lap with row equal to current node.

Definition at line 422 of file ImplicitVectorOperators_fwd.hpp.

◆ lap() [2/2]

template<class shape_storage_type , class matrix_type , class rhs_type >

BasicOp<Lap<dim> > mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::lap	(	int	node,
		int	row
	)

inline

Sets implicit equation that Laplacian of a vector field is equal to some other value.

The code alpha*op.lap(node) = v fills the node-th, node+N-th and node+2*N-th row of matrix M these rows of equation \( M u = v \) are a good approximation of the equation

\[ \alpha \nabla^2\vec{u}(p) = \vec{v}(p), \]

where \(p\) is the node-th point.

Note: Above example was made for dim = 3, other dimensions are analogous.

Parameters

node	Index of a node from 0 to `N` for which to write the equation for. This means only `node`-th`, node+N`-th and `node+2*N`-th rows of the matrix are changed. User must make sure that the matrix is large enough.
row	Write equation in this specific row. Row with index `node` is chosen by default.

Returns: A class representing this operation for lazy evaluation purposes.

Example:

    // Solve problem for F unknown vector field:
    // (Laplace F)(x, y) = [4 (1 - 2 pi^2 y^2) sin(2 pi x); 4 (1 - x - x y) e^(2 y)]
    // on [0, 1] x [0, 1] with BC:
    // F(x, 0) = [0, 1]
    // F(x, 1) = [2 sin(2 pi x), (1 - x) e^2]
    // F(0, y) = [0, e^(2 y)]
    // F(1, y) = [0, (1 - y) e^(2 y)]
    // with solution:
    // F(x, y) = [2 sin(2 pi x) y^2,  (1 - x y) e^(2 y)]
    std::function<Vec2d(Vec2d)> analytical = [] (const Vec2d& p) {
        return Vec2d({2*std::sin(2*PI*p[0])*p[1]*p[1],  (1 - p[0]*p[1]) * std::exp(2*p[1])});
    };
    // Prepare domain
    BoxShape<Vec2d> box(0., 1.);
    DomainDiscretization<Vec2d> domain = box.discretizeWithStep(0.025);
    int N = domain.size();
    int support_size = 9;
    domain.findSupport(FindClosest(support_size));
    // Prepare operators and matrix
    WLS<Monomials<Vec2d>, GaussianWeight<Vec2d>> wls(2, std::sqrt(15.0 / N));
    auto storage = domain.computeShapes(wls);
    Eigen::SparseMatrix<double> M(2*N, 2*N);
    M.reserve(storage.supportSizes()+storage.supportSizes());
    Eigen::VectorXd rhs(2*N); rhs.setZero();
    auto op = storage.implicitVectorOperators(M, rhs);
    // Set equation on interior
    for (int i : domain.interior()) {
        double x = domain.pos(i, 0), y = domain.pos(i, 1);
        op.lap(i) = {4 * (1 - 2 * PI*PI * y*y) * std::sin(2 * PI * x),
                     4 * (1 - x - x * y) * std::exp(2 * y)};
    }
    // Set boundary conditions
    for (int i : domain.types() == -1) {
        double y = domain.pos(i, 1);
        op.value(i) = {0, std::exp(2*y)};
    }
    for (int i : domain.types() == -2) {
        double y = domain.pos(i, 1);
        op.value(i) = {0, (1-y) * std::exp(2*y)};
    }
    for (int i : domain.types() == -3) {
        op.value(i) = {0, 1};
    }
    for (int i : domain.types() == -4) {
        double x = domain.pos(i, 0);
        op.value(i) = {2 * std::sin(2 * PI * x), (1 - x) * std::exp(2)};
    }
    // Sparse solve
    M.makeCompressed();
 
    Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
    solver.compute(M);
    Eigen::VectorXd sol = solver.solve(rhs);
    for (int i = 0; i < N; ++i) {
        Vec2d correct = analytical(domain.pos(i));
        EXPECT_NEAR(sol[i], correct[0], 1e-3);
        EXPECT_NEAR(sol[N+i], correct[1], 1e-3);
    }

Definition at line 420 of file ImplicitVectorOperators_fwd.hpp.

◆ neumann() [1/2]

template<class shape_storage_type , class matrix_type , class rhs_type >

GradOp mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::neumann	(	int	node,
		const domain_vector_t &	n
	)

inline

Same as neumann with row equal to current node.

Definition at line 511 of file ImplicitVectorOperators_fwd.hpp.

◆ neumann() [2/2]

template<class shape_storage_type , class matrix_type , class rhs_type >

GradOp mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::neumann	(	int	node,
		const domain_vector_t &	unit_normal,
		int	row
	)

inline

Sets neumann boundary conditions in node node.

The code alpha*op.neumann(node, normal) = v fills the node-th, node+N-th and node+2*N-th row of matrix M these rows of equation \( M u = v \) are a good approximation of the equation

\[ \alpha \dpar{\vec{u}}{\vec{n}}(p) = \vec{v}(p) \]

where \(p\) is the node-th point and \(\vec{n}\) is the unit_normal.

This is the same as using grad(), but has additional semantic meaning of setting the boundary conditions.

Returns: A class representing this operation for lazy evaluation purposes.

See also: grad

Definition at line 505 of file ImplicitVectorOperators_fwd.hpp.

◆ setColOffset()

template<class shape_storage_type , class matrix_type , class rhs_type >

void mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::setColOffset ( int col_offset )

inline

Sets col offset for given matrix, treating it as is the first column had index col_offset.

Definition at line 355 of file ImplicitVectorOperators_fwd.hpp.

◆ setProblem()

template<class shape_storage_type , class matrix_type , class rhs_type >

void mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::setProblem	(	matrix_t &	M,
		rhs_t &	rhs,
		int	row_offset = `0`,
		int	col_offset = `0`
	)

inline

Sets current matrix and right hand side.

Definition at line 343 of file ImplicitVectorOperators_fwd.hpp.

◆ setRowOffset()

template<class shape_storage_type , class matrix_type , class rhs_type >

void mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::setRowOffset ( int row_offset )

inline

Sets row offset for given matrix, treating it as is the first row had index row_offset.

Definition at line 350 of file ImplicitVectorOperators_fwd.hpp.

◆ traction() [1/2]

template<class shape_storage_type , class matrix_type , class rhs_type >

TractionOp mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::traction	(	int	node,
		scalar_t	lam,
		scalar_t	mu,
		const domain_vector_t &	n
	)

inline

Same as traction with row equal to current node.

Definition at line 539 of file ImplicitVectorOperators_fwd.hpp.

◆ traction() [2/2]

template<class shape_storage_type , class matrix_type , class rhs_type >

TractionOp mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::traction	(	int	node,
		scalar_t	lam,
		scalar_t	mu,
		const domain_vector_t &	unit_normal,
		int	row
	)

inline

Sets traction boundary conditions in node node.

The code alpha*op.traction(node, lam, mu, normal) = t fills the node-th, node+N-th and node+2*N-th row of matrix M these rows of equation \( M u = t \) are a good approximation of the equation

\[ \alpha (\sigma(\vec{u}) \vec{n}) (p) = \vec{t}(p) \]

where \(p\) is the node-th point, \(\sigma\) is the stress tensor and \(\vec{n}\) is the unit_normal.

Parameters

node	Index of a node from 0 to `N` for which to write the equation for. This means only `node`-th`, node+N`-th and `node+2*N`-th rows of the matrix are changed. User must make sure that the matrix is large enough.
row	Write equation in this specific row. Row with index `node` is chosen by default.
lam	The first Lame coefficient.
mu	The second Lame coefficient.
unit_normal	The outside unit normal on the surface for which to set the traction.

Returns: A class representing this operation for lazy evaluation purposes.

Definition at line 532 of file ImplicitVectorOperators_fwd.hpp.

◆ value() [1/2]

template<class shape_storage_type , class matrix_type , class rhs_type >

ValueOp mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::value ( int node )

inline

Same as value with row equal to current node.

Definition at line 398 of file ImplicitVectorOperators_fwd.hpp.

◆ value() [2/2]

template<class shape_storage_type , class matrix_type , class rhs_type >

ValueOp mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::value	(	int	node,
		int	row
	)

inline

Sets implicit equation that value of a vector field is equal to some other value.

The code alpha*op.value(node) = v fills the node-th, node+N-th and node+2*N-th row of matrix M these rows of equation \( M u = v \) are a good approximation of the equation

\[ \alpha \vec{u}(p) = \vec{v}(p), \]

where \(p\) is the node-th point.

Note: Above example was made for dim = 3, other dimensions are analogous.

Parameters

node	Index of a node from 0 to `N` for which to write the equation for. This means only `node`-th`, node+N`-th and `node+2*N`-th rows of the matrix are changed. User must make sure that the matrix is large enough.
row	Write equation in this specific row. Row with index `node` is chosen by default.

Returns: A class representing this operation for lazy evaluation purposes.

Definition at line 396 of file ImplicitVectorOperators_fwd.hpp.

Friends And Related Function Documentation

◆ operator<<

template<class shape_storage_type , class matrix_type , class rhs_type >

template<typename S , typename M , typename R >

std::ostream& operator<<	(	std::ostream &	os,
		const ImplicitVectorOperators< S, M, R > &	op
	)

friend

Output basic info about given operators.

Definition at line 43 of file ImplicitVectorOperators.hpp.

Member Data Documentation

◆ col_offset

template<class shape_storage_type , class matrix_type , class rhs_type >

int mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::col_offset

private

Column offset to be used when accessing matrix coefficients.

Definition at line 65 of file ImplicitVectorOperators_fwd.hpp.

◆ M

template<class shape_storage_type , class matrix_type , class rhs_type >

matrix_t* mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::M

private

Pointer to problem matrix.

Definition at line 62 of file ImplicitVectorOperators_fwd.hpp.

◆ rhs

template<class shape_storage_type , class matrix_type , class rhs_type >

rhs_t* mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::rhs

private

Pointer to right hand side.

Definition at line 63 of file ImplicitVectorOperators_fwd.hpp.

◆ row_offset

template<class shape_storage_type , class matrix_type , class rhs_type >

int mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::row_offset

private

Row offset to be used when accessing matrix or rhs coefficients.

Definition at line 64 of file ImplicitVectorOperators_fwd.hpp.

◆ ss

template<class shape_storage_type , class matrix_type , class rhs_type >

const shape_storage_t* mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >::ss

private

Shape storage, but name is shortened for readability.

Definition at line 61 of file ImplicitVectorOperators_fwd.hpp.

The documentation for this class was generated from the following files:

include/medusa/bits/operators/ImplicitVectorOperators_fwd.hpp
include/medusa/bits/operators/ImplicitVectorOperators.hpp

Detailed Description

template<class shape_storage_type, class matrix_type, class rhs_type> class mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >

Public Member Functions

Public Types

Classes

Friends

Private Member Functions

Private Attributes

Member Enumeration Documentation

◆ anonymous enum

Constructor & Destructor Documentation

◆ ImplicitVectorOperators() [1/3]

◆ ImplicitVectorOperators() [2/3]

◆ ImplicitVectorOperators() [3/3]

Member Function Documentation

◆ addToM()

◆ addToRhs()

◆ apply() [1/3]

◆ apply() [2/3]

◆ apply() [3/3]

◆ eq()

◆ grad() [1/2]

◆ grad() [2/2]

◆ graddiv() [1/2]

◆ graddiv() [2/2]

◆ hasMatrix()

◆ hasRhs()

◆ hasShapes()

◆ lap() [1/2]

◆ lap() [2/2]

◆ neumann() [1/2]

◆ neumann() [2/2]

◆ setColOffset()

◆ setProblem()

◆ setRowOffset()

◆ traction() [1/2]

◆ traction() [2/2]

◆ value() [1/2]

◆ value() [2/2]

Friends And Related Function Documentation

◆ operator<<

Member Data Documentation

◆ col_offset

◆ M

◆ rhs

◆ row_offset

◆ ss

template<class shape_storage_type, class matrix_type, class rhs_type>
class mm::ImplicitVectorOperators< shape_storage_type, matrix_type, rhs_type >