(Bi)linear forms#

Bilinear forms#

Bilinear forms are managed using:

Figure made with TikZ

Figure made with TikZ

The BasicBilinearForm class#

The BasicBilinearForm class is an abstract class having only the unknowns as attribute and their associated accessors:

class BasicBilinearForm
{
  protected:
    const Unknown* u_p; // pointer to unknown u
    const Unknown* v_p; // pointer to test function v
  public:
    const Unknown& up() const;
    const Unknown& vp() const;
  ...

The u letter always refers to the left unknown of the bilinear form and the v letter to the right unknown (test function) of the bilinear form. In matrix representation u stands for columns and v for rows.

The child classes inheriting from BasicBilinearForm have to provide the following member functions:

virtual ~BasicBilinearForm() {}
virtual BasicBilinearForm* clone() const = 0;
virtual LinearFormType type() const = 0;
virtual ValueType valueType() const = 0;
virtual void print(std::ostream&) const = 0;

The clone function is used to construct a copy of the child objects from parent.

The IntgBilinearForm class#

The IntgBilinearForm class handles linear forms defined by a (single) integral over a geometric domain:

\[\int_{D} \mathcal{L}_1(u) op \mathcal{L}_2(v)\]

where \(D\) is a geometric domain (GeomDomain class) and \({\cal L}_1\), \({\cal L}_2\) are linear operators (OperatorOnUnknown class) acting on unknowns \(u\) and \(v\) (Unknown class) and \(op\) is an algebraic product operator. In case of mesh domain, it handles also the quadrature rules, one rule by shape of elements.

This class is defined as follows:

class IntgBilinearForm : public BasicBilinearForm
{
  protected:
    const GeomDomain* domain_p; // geometric domain of the integral (pointer)
    const OperatorOnUnknown* opu_p; // operator on unknown u (pointer)
    const OperatorOnUnknown* opv_p; // operator on unknown v (pointer)
    AlgebraicOperator aop_;         // algebraic operation between operators
    map<Number, Quadrature*> quadratures_; // pointers to quadrature rules
    ...

Note that the quadratures_ map may be empty. This is the case when the domain is not a mesh domain. Generally, the elements of a mesh domain have the same shape, but it may not. It is the reason why a map is used to store the quadratures.

It has one basic constructor, some public accessors and utilities:

IntgBilinearForm(const GeomDomain&,const OperatorOnUnknown&, AlgebraicOperator,
                const OperatorOnUnknown&, QuadRule = _defaultRule, Number = 0)
const OperatorOnUnknown* opu() const;
const OperatorOnUnknown* opv() const;
const GeomDomain* domain() const;
virtual BasicBilinearForm* clone() const;
virtual LinearFormType type() const;
virtual ValueType valueType() const;
void setQuadrature(QuadRule, Number);
virtual void print(std::ostream&) const;

The setQuadrature function set the quadrature pointers from a QuadRule (enumeration of type of quadrature rule) and a quadrature number (degree or order of quadrature rule). For a full description of quadrature rules see section The Quadrature and the QuadratureRule classes. When QuadRule is set to _defaultRule in constructor, the ‘best’ rule is chosen regarding the shape element and the degree of integrand (see Quadrature::bestQuadRule function).

The DoubleIntgBilinearForm class#

The DoubleIntgBilinearForm class handles linear forms defined by a double integral over a product of geometric domains:

\[\int_{D_x}\int_{D_y} {\cal L}(u)\]

where \(D_x\), \(D_y\) are geometric domain (GeomDomain class) and \({\cal L}\) a linear operator (OperatorOnUnknown class) acting on unknown \(u\) (Unknownn class). Except for the fact that there are two geometric domains, this class is very similar to the IntgBilinearForm class:

class DoubleIntgBilinearForm : public BasicBilinearForm
{
  protected:
    const GeomDomain* domainx_p;    // first geometric domain (say x variable)
    const GeomDomain* domainy_p;    // second geometric domain (say y variable)
    const OperatorOnUnknown* opu_p; // operator on unknown u (pointer)
    const OperatorOnUnknown* opv_p; // operator on unknown v (pointer)
    ...

It provides a basic constructor and some public accessors:

DoubleIntgBilinearForm(const GeomDomain&, const GeomDomain&, const OperatorOnUnknown&, AlgebraicOperator, const OperatorOnUnknown&);
const OperatorOnUnknown* opu() const;
const OperatorOnUnknown* opv() const;
const GeomDomain* domainx() const;
const GeomDomain* domainy() const;
virtual BasicBilinearForm* clone() const;
virtual LinearFormType type() const;
virtual ValueType valueType() const;
virtual void print(std::ostream&) const;

Warning

Be cautious in the definition of double integral. The syntax:

BilinearForm a=intg(Sigma, Gamma ,u*G*v);

has to be interpreted as

\[\int_{\Sigma_x}\int_{\Gamma_y}u(y)G(x,y)v(x)dy\,dx\]

where \(u\) (the right unknown) stands for the unknown, \(v\) (the left unknown) being the test function. In matrix representation, it means that \(u\) stands for columns and \(v\) for rows.

The UserBilinearForm class#

The UserBilinearForm class allows users to define general bilinear form. It is based on a BlfFunction providing the computation of elementary matrices. This function has always the following signature:

void blfun(BlfDataComputation& blfd);

where the BlfDataComputation class encapsulates some useful data updated by computation algorithms:

class BlfDataComputation
{
  public:
    const Element* elt_u, *elt_v;     // Element pointers (FEM or BEM)
    const Element* elt_u2,*elt_v2;    // additional Element pointers (DG)
    const GeomElement* sidelt;        // side element when DG
    vector<Matrix<Real>> matels;      // real elementary matrices
    vector<Matrix<Complex>> cmatels;  // complex elementary matrices
    ...

Regarding their own business, users have to fill either real elementary matrices (matels) or complex elementary matrices (cmatels) from element data (elt_u, elt_v, …).

Caution

Up to now, only FEM and DG computation algorithms manage UserBilinearForm.

So, the UserBilinearForm class handles the following data:

class UserBilinearForm : public BasicBilinearForm
{
  public:
    BlfFunction blfun_;                    // pointer to an external blf function (0 by default)
    const IntegrationMethod* intgMethod_p; // pointer to an integration method
    bool requireInvJacobian_;              // requiring jacobian (default=false)
    bool requireNormal_;                   // requiring normal vector (default=false)

This class provides two general constructors and a copy constructor:

UserBilinearForm(const GeomDomain& dom, const Unknown& u, const Unknown& v, BlfFunction blf, ComputationType ct,SymType st, bool rij, bool rno, const IntegrationMethod& im);
UserBilinearForm(const GeomDomain& domv, const GeomDomain& domu, const Unknown& u, const Unknown& v, BlfFunction blf, ComputationType ct, SymType st, bool rij, bool rno, const IntegrationMethod& im);
UserBilinearForm(const UserBilinearForm&);
BasicBilinearForm* clone() const; // clone of the UserBilinearForm
~UserBilinearForm();

In practice, it is advised to use external pseudo-constructors with some default arguments:

BilinearForm userBlf(const GeomDomain& g, const Unknown& u, const Unknown& v, BlfFunction blf, ComputationType ct, SymType=_noSymmetry, bool reqIJ=false, bool reqN=false, const IntegrationMethod& im=QuadratureIM(_GaussLegendreRule,3));
BilinearForm userBlf(const GeomDomain& g, const GeomDomain&, const Unknown& u, const Unknown& v, BlfFunction blf, ComputationType ct, SymType=_noSymmetry, bool reqIJ=false, bool reqN=false, const IntegrationMethod& im=QuadratureIM(_GaussLegendreRule,3));

Besides some useful tools are provided:

LinearFormType type() const;                 // return the type of the linear form (_userLf)
const IntegrationMethod* intgMethod() const; // pointer to the integration method object
ValueType valueType() const;                 // return the value type
void requireInvJacobian(bool tf=true);       // set on/off the requireInvJacobian flag
void requireNormal(bool tf=true);            // set on/off the requireNormal flag
SymType setSymType() const;                  // set the symmetry property (done by analysis)
void print(std::ostream&) const;             // print utility
void print(PrintStream& os) const;           // print utility
String asString() const;                     // interpret as string for print purpose

The next simple example shows how to deal with UserBilinearForm. It deals with the \(\int_\Omega\nabla u.\nabla v\) bilinear form that can be handled by standard approach:

// compute elementary matrix grad.grad in 2D-P1
void blfGradGrad(BlfDataComputation& blfd)
{
  if (blfd.matel().size() == 0) blfd.matel() = Matrix<Real>(3,3);
  const Element* elt=blfd.elt_u; // get element concerned by
  if (elt == 0) return;          // no computation
  GeomMapData& mapdata = *elt->geomElt_p->meshElement()->geomMapData_p; // get geometric data
  Matrix<Real> C=(0.5*mapdata.differentialElement) * mapdata.inverseJacobianMatrix * tran(mapdata.inverseJacobianMatrix);
  Matrix<Real> G(2,3,0.); G(1,1)=1; G(2,2)=1; G(1,3)=-1; G(2,3)=-1; // gradient in ref element
  blfd.matel()=tran(G)*C*G; // elementary matrix
}
...
Mesh ms(Rectangle(_xmin=0., _xmax=1., _ymin=0., _ymax=1., _nnodes=10, _domain_name="Omega"),
        _shape=_triangle, _order=1, _generator=_structured, _split_direction=_alternate);
Domain omega=ms.domain("Omega");
Space H(_domain=omega, _interpolation=P1, _name="H");
Unknown p(H, _name="p"); TestFunction t(p, _name="t");
BilinearForm ublf = userBlf(omega,p,t,blfGradGrad,_FEComputation,_symmetric, true, false); // define user blf
TermMatrix Ku(ublf, _name="Ku"); // use it as usual

See also

See the unit\_DG.cpp file for a more complex example which deals with the discontinuous Galerkin term:

\[\int_\Gamma \big\{\nabla_hu_h.n\big\}\big[v\big]\]

The SuBilinearForm class#

BasicBilinearForm objects may be linearly combined to produce a linear combination of BasicBilinearForm objects stored as a list of pair of BasicBilinearForm object and a complex scalar in the SuBilinearForm class.

typedef std::pair<BasicBilinearForm*,Complex> blfPair;
class SuBilinearForm
{
  protected:
    std::vector<blfPair> blfs_; //list of pairs of basic bilinear form and coefficient
    ...

All BasicBilinearForm objects must have the same pair of unknowns! It is the reason why this class has no pointer to unknowns; it refers to the first basic bilinear form to get its unknowns.

Because BasicBilinearForm objects are copied for safety reason, this class provides default and basic constructors but also a copy constructor, a destructor and the overload assignment operator:

SuBilinearForm() {}
SuBilinearForm(const SuBilinearForm&);
~SuBilinearForm();
SuBilinearForm& operator=(const SuBilinearForm&);

It provides few accessors (const and non const):

Number size() const;
std::vector<blfPair>& blfs();
const std::vector<blfPair>& blfs() const;
blfPair& operator()(Number n);
const blfPair& operator()(Number n) const;
const Unknown* up() const;
const Unknown* vp() const;
const Space* uSpace() const;
const Space* vSpace() const;
LinearFormType type() const;
ValueType valueType() const;

It is possible to perform linear combination of linear combinations using the following overloaded operators:

SuBilinearForm& SuBilinearForm::operator +=(const SuBilinearForm&);
SuBilinearForm& SuBilinearForm::operator -=(const SuBilinearForm&);
SuBilinearForm& SuBilinearForm::operator *=(const Complex&);
SuBilinearForm& SuBilinearForm::operator /=(const Complex&);
SuBilinearForm operator-(const SuBilinearForm&);
SuBilinearForm operator+(const SuBilinearForm&, const SuBilinearForm&);
SuBilinearForm operator-(const SuBilinearForm&, const SuBilinearForm&);
SuBilinearForm operator*(const Complex&, const SuBilinearForm&);
SuBilinearForm operator*(const SuBilinearForm&, const Complex&);
SuBilinearForm operator/(const SuBilinearForm&, const Complex&);
bool SuBilinearForm::checkConsistancy(const SuBilinearForm&) const;

The member function checkConsistancy performs a test to ensure that pair of unknowns is always the same.

Finally, there are some print facilities:

void SuBilinearForm::print(std::ostream&) const;
std::ostream& operator<<(std::ostream&,const SuBilinearForm&);

The BilinearForm class#

The BilinearForm class is the end user class dealing with general bilinear form, either a single unknown bilinear form (a SuBilinearForm object) or a multiple unknown bilinear form (a list of SuBilinearForm objects). In this class, a single unknown bilinear form is a multiple unknown bilinear form with one pair of unknowns! The list of SuBilinearForm objects is stored in a map of SuBilinearForm, indexed by a pair of pointers to unknown:

typedef std::pair<const Unknown *,const Unknown *> uvPair;

class BilinearForm
{
  protected:
    std::map<uvPair,SuBilinearForm> mlcblf_; // list of linear combinations of basic forms
    ...

To manage the map, the following aliases are defined:

typedef std::map<uvPair,SuBilinearForm>::iterator it_mublc;
typedef std::map<uvPair,SuBilinearForm>::const_iterator cit_mublc;

Caution

When the SuBilinearForm unknowns own an unknown component, the SuBilinearForm object is attached to its parent item! In other words, component unknowns are not indexed in the map.

This class provides only one constructor from a linear combination of forms:

BilinearForm(const SuBilinearForm&);

and proposes some accessors and facilities:

bool singleUnknown() const;
bool isEmpty() const;
SuBilinearForm& operator[](const uvPair&);             // protected
const SuBilinearForm& operator[](const uvPair&) const; // protected
const SuBilinearForm& first() const;
BilinearForm operator()(const Unknown&, const Unknown&) const;
BasicBilinearForm& operator()(const Unknown&, const Unknown&,Number);
const BasicBilinearForm& operator()(const Unknown&, const Unknown&, Number) const;

Besides, there is the overloaded external end-user’s function intg that constructs BilinearForm object:

BilinearForm intg(const GeomDomain&, const OperatorOnUnknown&, AlgebraicOperator, const OperatorOnUnknown&);
BilinearForm intg(const GeomDomain&, const OperatorOnUnknowns&);
BilinearForm intg(const GeomDomain&, const LcOperatorOnUnknowns&);
BilinearForm intg(const GeomDomain&, const GeomDomain&, const OperatorOnUnknown&, AlgebraicOperator, const OperatorOnUnknown&);
BilinearForm intg(const GeomDomain&, const GeomDomain&, const OperatorOnUnknowns&);

OperatorUnknowns is a simple class used to store the two operators and the algebraic operation.

In order to construct any linear forms, the algebraic operators (+=, -=, *=, /=, +, -, *, /) are oveloaded for different objects:

BilinearForm& BilinearForm::operator+=(const BilinearForm&);
BilinearForm& BilinearForm::operator-=(const BilinearForm&);
BilinearForm& BilinearForm::operator*=(const Complex&);
BilinearForm& BilinearForm::operator/=(const Complex&);
BilinearForm operator+(const BilinearForm&, const BilinearForm&);
BilinearForm operator-(const BilinearForm&, const BilinearForm&);
BilinearForm operator*(const Complex&, const BilinearForm&);
BilinearForm operator*(const BilinearForm&, const Complex&);
BilinearForm operator/(const BilinearForm&, const Complex&);

Finally, the class provides usual print facilities:

void BilinearForm::print(std::ostream&) const;
std::ostream& operator<<(std::ostream&, const BilinearForm&);

Example#

To end we show some characteristic examples. In these examples, u, v are scalar unknowns, p, q are vector unknowns, f either a scalar or a scalar function, h either a vector or a vector function, g a scalar or a scalar kernel and Omega, Gamma, Sigma some geometric domains.

BilinearForm a1 = intg(Omega, grad(u)|grad(v));
a1 -= k2*intg(Omega, f*u*v);
a1 = intg(Omega, grad(u)|grad(v)) - intg(Omega, f*u*v); // same result
BilinearForm a2 = intg(Omega, div(p)*q) + eps*intg(Omega, p*q);
BilinearForm a3 = a1 + a2;
a3 = intg(Omega, grad(u)|grad(v)) + intg(Omega, div(p)*q) - intg(Omega, f*u*v) + eps*intg(Gamma, p*q); // same result

Linear forms#

Linear forms are managed using:

Figure made with TikZ

Figure made with TikZ

The BasicLinearForm class#

The BasicLinearForm class is an abstract class having only an unknown as attribute and its associated accessor:

class BasicLinearForm
{
  protected:
    const Unknown* u_p;       // pointer to unknown
  public:
    const Unknown& up() const // return the unknown

The child classes inheriting from BasicLinearForm have to provide the following member functions:

virtual ~BasicLinearForm() {}                // virtual destructor
virtual BasicLinearForm* clone() const = 0;  // clone of the linear form
virtual LinearFormType type() const = 0;     // return the type of the linear form
virtual ValueType valueType() const = 0;     // return the value type of the linear form
virtual void print(std::ostream&) const = 0; // print utility

The virtual clone function is used to construct a copy of the child objects from parent.

The IntgLinearForm class#

The IntgLinearForm class handles linear forms defined by a (single) integral over a geometric domain:

\[\int_{D} {\cal L}(u)\]

where \(D\) is a geometric domain (GeomDomain class) and \({\cal L}\) a linear operator (OperatorOnUnknown class) acting on unknown \(u\) (Unknown class).

This class is defined as follows:

class IntgLinearForm : public BasicLinearForm
{
  protected:
    const GeomDomain* domain_p; // geometric domain of the integral (pointer)
    const OperatorOnUnknown* opu_p; // operator on unknown (pointer)

It has some public constructors and some public accessors:

IntgLinearForm(const GeomDomain&, const OperatorOnUnknown&);
const OperatorOnUnknown* opu() const;
const GeomDomain* domain() const;
virtual BasicLinearForm* clone() const;
virtual LinearFormType type() const
virtual ValueType valueType() const
virtual void print(std::ostream&) const;

The DoubleIntgLinearForm class#

The DoubleIntgLinearForm class handles linear forms defined by a double integral over a product of geometric domains:

\[\int_{D_x}\int_{D_y} \mathcal{L}(u)\]

where \(D_x\), \(D_y\) are geometric domain (GeomDomain class) and \(\mathcal{L}\) a linear operator (OperatorOnUnknown class) acting on unknown \(u\) (Unknown class).

Except for the fact that there are two geometric domains, this class is very similar to the IntgLinearForm class:

class DoubleIntgLinearForm : public BasicLinearForm
{
  protected:
  const OperatorOnUnknown* opu_p; // operator on unknown (pointer)
  const GeomDomain* domainx_p; // first geometric domain (say x variable)
  const GeomDomain* domainy_p; // second geometric domain (say y variable)

It provides some public constructors and some public accessors:

DoubleIntgLinearForm(const GeomDomain&, const GeomDomain&,const OperatorOnUnknown&);
const OperatorOnUnknown* opu() const;
const GeomDomain* domainx() const;
const GeomDomain* domainy() const;
virtual BasicLinearForm* clone() const;
virtual LinearFormType type() const;
virtual ValueType valueType() const;
virtual void print(std::ostream&) const;

The SuLinearForm class#

BasicLinearForm objects may be linearly combined to produce a linear combination of BasicLinearForm objects stored as a list of pair of BasicLinearForm object and a complex scalar in the SuLinearForm class.

typedef std::pair<BasicLinearForm*,Complex> lfPair;
class SuLinearForm
{
  protected:
    std::vector<lfPair> lfs_; //list of pairs of basic linear form and coefficient
    ...

All BasicLinearForm objects must have the same unknown! It is the reason why this class has no pointer to an unknown; it refers to the first basic linear form to get its unknown.

Because BasicLinearForm objects are copied for safety reason, this class provides default and basic constructors but also a copy constructor, a destructor and the overload assignment operator:

SuLinearForm() {}
SuLinearForm(const SuLinearForm&);
~SuLinearForm();
SuLinearForm& operator=(const SuLinearForm&);

It provides few accessors (const and non-const):

Number size() const;
std::vector<lfPair>& lfs();
const std::vector<lfPair>& lfs() const;
lfPair& operator()(Number n);
const lfPair& operator()(Number n) const;
const Unknown* unknown() const;
const Space* space() const;
LinearFormType type() const;
ValueType valueType() const;

It is possible to perform linear combination of linear combinations using the following overloaded operators:

SuLinearForm& SuLinearForm::operator +=(const SuLinearForm&);
SuLinearForm& SuLinearForm::operator -=(const SuLinearForm&);
SuLinearForm& SuLinearForm::operator *=(const Complex&);
SuLinearForm& SuLinearForm::operator /=(const Complex&);
SuLinearForm operator-(const SuLinearForm&);
SuLinearForm operator+(const SuLinearForm&, const SuLinearForm&);
SuLinearForm operator-(const SuLinearForm&, const SuLinearForm&);
SuLinearForm operator*(const Complex&, const SuLinearForm&);
SuLinearForm operator*(const SuLinearForm&, const Complex&);
SuLinearForm operator/(const SuLinearForm&, const Complex&);
bool SuLinearForm::checkConsistancy(const SuLinearForm&) const;

The member function checkConsistancy performs a test to ensure that the unknown is always the same.

Finally, there are some print facilities:

void SuLinearForm::print(std::ostream&) const;
std::ostream& operator<<(std::ostream&,const SuLinearForm&);

The LinearForm class#

The LinearForm class is the end-user class dealing with general linear form, either a single unknown linear form (a SuLinearForm object) or a multiple unknown linear form (a list of SuLinearForm objects). In this class, a single unknown linear form is a multiple unknown linear form with one unknown! The list of SuLinearForm objects is stored in a map of SuLinearForm, indexed by the pointer to SuLinearForm unknown:

class LinearForm
{
  protected:
    std::map<const Unknown *,SuLinearForm> mlclf_; // list of linear combinations of basic forms
    ...

To manage the map, the following aliases are defined:

typedef std::map<const Unknown *,SuLinearForm>::iterator it_mulc;
typedef std::map<const Unknown *,SuLinearForm>::const_iterator cit_mulc;

Caution

When the SuLinearForm unknown is an unknown component, the SuLinearForm object is attached to its parent item! In other words, component unknowns is not indexed in the map.

This class provides only one constructor from a linear combination of forms:

LinearForm(const SuLinearForm&);

and offers some accessors and facilities:

bool isEmpty() const;
bool singleUnknown() const;
SuLinearForm& operator[](const Unknown*);             // protected
const SuLinearForm& operator[](const Unknown*) const; // protected
const SuLinearForm& first() const;
LinearForm operator()(const Unknown&) const;
BasicLinearForm& operator()(const Unknown&, Number);
const BasicLinearForm& operator()(const Unknown&, Number) const;

Besides, there is the overloaded external end-user’s function (intg) that constructs LinearForm object:

LinearForm intg(const GeomDomain&, const LcOperatorOnUnknown&);
LinearForm intg(const GeomDomain&, const OperatorOnUnknown&);
LinearForm intg(const GeomDomain&, const Unknown&);
LinearForm intg(const GeomDomain&, const GeomDomain&, const OperatorOnUnknown&);
LinearForm intg(const GeomDomain&, const GeomDomain&, const Unknown&);

In order to construct any linear forms, the algebraic operators (+=, -=, *=, /=, +, -, *, /) are overloaded for different objects:

LinearForm& LinearForm::operator +=(const LinearForm&);
LinearForm& LinearForm::operator -=(const LinearForm&);
LinearForm& LinearForm::operator *=(const Complex&);
LinearForm& LinearForm::operator /=(const Complex&);
LinearForm operator+(const LinearForm&, const LinearForm&);
LinearForm operator-(const LinearForm&, const LinearForm&);
LinearForm operator*(const Complex&, const LinearForm&);
LinearForm operator*(const LinearForm&, const Complex&);
LinearForm operator/(const LinearForm&, const Complex&);

Finally, the class provides usual print facilities:

void LinearForm::print(std::ostream&) const;
std::ostream& operator<<(std::ostream&, const LinearForm&);

Example#

To end we show some characteristic examples. In these examples, u is a scalar unknown, w is a vector unknown, f either a scalar or a scalar function, h either a vector or a vector function, g a scalar or a scalar kernel and Omega,:xlifepp:Gamma Sigma some geometric domains.

LinearForm l1=intg(Omega, f*u);
l1 += intg(Omega, h|grad(u));
l1 = intg(Omega, f*u) + intg(Omega, h|grad(u)); // equivalent
l1 = l1 - 2*intg(Gamma, Sigma, g*u);
LinearForm l2 = intg(Omega, f*div(w)) + intg(Gamma, h|(w^n));
LinearForm l3 = l1 + l2; // multiple unknowns
l3 = intg(Omega, f*u) + intg(Omega, f*div(w)) + intg(Omega, h|grad(u)) + 2*intg(Gamma, Sigma, g*u) + intg(Gamma, h|(w^n)); // same result