diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 13a6d08a..979e95ef 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.11.1","generation_timestamp":"2024-10-17T13:11:13","documenter_version":"1.7.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.11.1","generation_timestamp":"2024-10-22T11:52:01","documenter_version":"1.7.0"}} \ No newline at end of file diff --git a/dev/assemble/index.html b/dev/assemble/index.html index 33689402..e7633aee 100644 --- a/dev/assemble/index.html +++ b/dev/assemble/index.html @@ -58,4 +58,4 @@ ) end -integrand(op::SingleLayerTrace, kernel, g, τ, f, σ) = f[1]*g[1]*kernel.green

Every kernel corresponds with a type. Kernels can potentially depend on a set of parameters; these appear as fields in the type. Here our Nitsche kernel depends on the wavenumber. In quaddata we precompute quadrature points for all geometric cells in the supports of test and trial elements. This is fairly sloppy: only one rule for test and trial integration is considered. A high accuracy implementation would typically compute points for both low quality and high quality quadrature rules.

Also quadrule is sloppy: we always select a DoubleQuadRule to perform the computation of interactions between local shape functions. No singularity extraction or other advanced technique is considered for nearby interactions. Clearly amateurs at work here!

BEAST provides a default implementation of an integration routine using double numerical quadrature. All that is required to tap into that implementation is a method overloading integrand. From the above formula it is clear what this method should look like.

That's it!

+integrand(op::SingleLayerTrace, kernel, g, τ, f, σ) = f[1]*g[1]*kernel.green

Every kernel corresponds with a type. Kernels can potentially depend on a set of parameters; these appear as fields in the type. Here our Nitsche kernel depends on the wavenumber. In quaddata we precompute quadrature points for all geometric cells in the supports of test and trial elements. This is fairly sloppy: only one rule for test and trial integration is considered. A high accuracy implementation would typically compute points for both low quality and high quality quadrature rules.

Also quadrule is sloppy: we always select a DoubleQuadRule to perform the computation of interactions between local shape functions. No singularity extraction or other advanced technique is considered for nearby interactions. Clearly amateurs at work here!

BEAST provides a default implementation of an integration routine using double numerical quadrature. All that is required to tap into that implementation is a method overloading integrand. From the above formula it is clear what this method should look like.

That's it!

diff --git a/dev/assets/postproc/index.html b/dev/assets/postproc/index.html index 7aa16b14..5f9cbb9d 100644 --- a/dev/assets/postproc/index.html +++ b/dev/assets/postproc/index.html @@ -1,4 +1,4 @@ Post-processing and Visualisation · BEAST Documentation
GitHub

Post-processing and Visualisation

Field computation

The main API for the computation of fields in a vector of points is the potential function. This function is capable of computing

\[F(x) = \int_{\Gamma} K(x,y) u(y) dy\]

with $u(y) = \Sigma_{i=1}^N u_i f_i(y)$ and $K(x,y)$ the integation kernel defining the type of potential. For example, the far field for a vector valued surface density is

\[F(x) = \int_{\Gamma} e^{ik \frac{x \cdot y}{|x|}} u(y) dy\]

For a far field potential such as this, the value only depends on the direction, not the magnitude, of $x$, as can be read off from the normalisation in the exponent.

The following script computes the far field along a semi-circle in the xz-plane.

Θ, Φ = range(0.0,stop=2π,length=100), 0.0
 dirs = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for θ in Θ for ϕ in Φ]
-farfield = potential(MWFarField3D(wavenumber=κ), dirs, u, X)
+farfield = potential(MWFarField3D(wavenumber=κ), dirs, u, X) diff --git a/dev/index.html b/dev/index.html index b487e085..292d946a 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,3 +1,3 @@ -BEAST.jl documentation · BEAST Documentation
GitHub

BEAST.jl documentation

BEAST provides a number of types modelling concepts and a number of algorithms for the efficient and simple implementation of boundary and finite element solvers. It provides full implementations of these concepts for the LU based solution of boundary integral equations for the Maxwell and Helmholtz systems.

Because Julia only compiles code at execution time, users of this library can hook into the code provided in this package at any level. In the extreme case it suffices to provide overwrites of the assemble functions. In that case, only the LU solution will be performed by the code here.

At the other end it suffices that users only supply integration kernels that act on the element-element interaction level. This package will manage all required steps for matrix assembly.

For the Helmholtz 2D and Maxwell 3D systems, complete implementations are supplied. These models will be discussed in detail to give a more concrete idea of the APIs provides and how to extend them.

Central to the solution of boundary integral equations is the assembly of the system matrix. The system matrix is fully determined by specifying a kernel G, a set of trial functions, and a set of test functions.

Basis

Sets of both trial and testing functions are implemented by models following the basis concept. The term basis is somewhat misleading as it is nowhere required nor enforced that these functions are linearly independent. Models implementing the Basis concept need to comply to the following semantics.

  • numfunctions(basis): number of functions in the Basis.
  • coordtype(basis): type of (the components of) the values taken on by the functions in the Basis.
  • scalartype(d): the scalar field underlying the vector space the basis functions take value in.
  • refspace(basis): returns the ReferenceSpace of local shape functions on which the Basis is built.
  • assemblydata(basis): assemblydata returns an iterable collection elements of geometric elements and a look table ad for use in assembly of interaction matrices. In particular, for an index element_idx into elements and an index local_shape_idx in basis of local shape functions refspace(basis), ad[element_idx, local_shape_idx] returns the iterable collection of (global_idx, weight) tuples such that the local shape function at local_shape_idx defined on the element at element_idx contributes to the basis function at global_idx with a weight of weight.
  • geometry(basis): returns an iterable collection of Elements. The order in which these Elements are encountered corresponds to the indices used in the assembly data structure.

Reference Space

The reference space concept defines an API for working with spaces of local shape functions. The main role of objects implementing this concept is to allow specialization of the functions that depend on the precise reference space used.

The functions that depend on the type and value of arguments modeling reference space are:

Kernel

A kernel is a fairly simple concept that mainly exists as part of the definition of a Discrete Operator. A kernel should obey the following semantics:

In many function definitions the kernel object is referenced by operator or something similar. This is a misleading name as an operator definition should always be accompanied by the domain and range space.

Discrete Operator

Informally speaking, a Discrete Operator is a concept that allows for the computation of an interaction matrix. It is a kernel together with a test and trial basis. A Discrete Operator can be passed to assemble and friends to compute its matrix representation.

A discrete operator is a triple (kernel, test_basis, trial_basis), where kernel is a Kernel, and test_basis and trial_basis are Bases. In addition, the following expressions should be implemented and behave according to the correct semantics:

In the context of fast methods such as the Fast Multipole Method other algorithms on Discrete Operators will typically be defined to compute matrix vector products. These algorithms do not explicitly compute and store the interaction matrix (this would lead to unacceptable computational and memory complexity).

BEAST.elements — Function

elements(geo)

Create an iterable collection of the elements stored in geo. The order in which this collection produces the elements determines the index used for lookup in the data structures returned by assemblydata and quaddata.

source
BEAST.scalartype — Function
scalartype(x)

The scalar field over which the values of a global or local basis function, or an operator are defined. This should always be a scalar type, even if the basis or operator takes on values in a vector or tensor space. This data type is used to determine the eltype of assembled discrete operators.

source
BEAST.assemblydata — Function
charts, admap, act_to_global = assemblydata(basis; onlyactives=true)

Given a basis this function returns a data structure containing the information required for matrix assemble, that is, the vector charts containing Simplex elements, a variable admap of type AssemblyData, and a mapping from indices of actively used simplices to global simplices.

When onlyactives is true, another layer of indices is introduced to filter out all cells of the mesh that are not in the union of the support of the basis functions (i.e., when the basis functions are defined only on a part of the mesh).

admap is, in essence, a three-dimensional array of named tuples, which, by wrapping it in the struct AssemblyData, allows the definition of iterators. The tuple consists of the two entries

admap[i,r,c].globalindex
-admap[i,r,c].coefficient

Here, c and r are indices in the iterable set of (active) simplices and the set of shape functions on each cell/simplex: r ranges from 1 to the number of shape functions on a cell/simplex, c ranges from 1 to the number of active simplices, and i ranges from 1 to the number of maximal number of basis functions, where any of the shape functions contributes to.

For example, for continuous piecewise linear lagrange functions (c0d1), each of the three shape functions on a triangle are associated with exactly one Lagrange function, and therefore i is limited to 1.

Note: When onlyactives=false, the indices c correspond to the position of the corresponding cell/simplex whilst iterating over geometry(basis). When onlyactives=true, then act_to_global(c) correspond to the position of the corresponding cell/simplex whilst iterating over geometry(basis).

For a triplet (i,r,c), globalindex is the index in the basis of the ith basis function that has a contribution from shape function r on (active) cell/simplex c. coefficient is the coefficient of that contribution in the linear combination defining that basis function in terms of shape function.

source
BEAST.geometry — Function
geometry(basis)

Returns an iterable collection of geometric elements on which the functions in basis are defined. The order the elements are encountered needs correspond to the element indices used in the data structure returned by assemblydata.

source
BEAST.refspace — Function
refspace(basis)

Returns the ReferenceSpace of local shape functions on which the basis is built.

source
BEAST.quaddata — Function
quaddata(operator, test_refspace, trial_refspace, test_elements, trial_elements)

Returns an object cashing data required for the computation of boundary element interactions. It is up to the client programmer to decide what (if any) data is cached. For double numberical quadrature, storing the integration points for example can significantly speed up matrix assembly.

  • operator is an integration kernel.
  • test_refspace and trial_refspace are reference space objects. quadata

is typically overloaded on the type of these local spaces of shape functions. (See the implementation in maxwell.jl for an example).

  • test_elements and trial_elements are iterable collections of the geometric

elements on which the finite element space are defined. These are provided to allow computation of the actual integrations points - as opposed to only their coordinates.

source
BEAST.quadrule — Function
quadrule(operator, test_refspace, trial_refspace, test_index, test_chart, trial_index, trial_chart, quad_data)

Based on the operator kernel and the test and trial elements, this function builds an object whose type and data fields specify the quadrature rule that needs to be used to accurately compute the interaction integrals. The quad_data object created by quaddata is passed to allow reuse of any precomputed data such as quadrature points and weights, geometric quantities, etc.

The type of the returned quadrature rule will help in deciding which method of momintegrals to dispatch to.

source
BEAST.momintegrals! — Function

momintegrals!(biop, tshs, bshs, tcell, bcell, interactions, strat)

Function for the computation of moment integrals using simple double quadrature.

source
+BEAST.jl documentation · BEAST Documentation
GitHub

BEAST.jl documentation

BEAST provides a number of types modelling concepts and a number of algorithms for the efficient and simple implementation of boundary and finite element solvers. It provides full implementations of these concepts for the LU based solution of boundary integral equations for the Maxwell and Helmholtz systems.

Because Julia only compiles code at execution time, users of this library can hook into the code provided in this package at any level. In the extreme case it suffices to provide overwrites of the assemble functions. In that case, only the LU solution will be performed by the code here.

At the other end it suffices that users only supply integration kernels that act on the element-element interaction level. This package will manage all required steps for matrix assembly.

For the Helmholtz 2D and Maxwell 3D systems, complete implementations are supplied. These models will be discussed in detail to give a more concrete idea of the APIs provides and how to extend them.

Central to the solution of boundary integral equations is the assembly of the system matrix. The system matrix is fully determined by specifying a kernel G, a set of trial functions, and a set of test functions.

Basis

Sets of both trial and testing functions are implemented by models following the basis concept. The term basis is somewhat misleading as it is nowhere required nor enforced that these functions are linearly independent. Models implementing the Basis concept need to comply to the following semantics.

  • numfunctions(basis): number of functions in the Basis.
  • coordtype(basis): type of (the components of) the values taken on by the functions in the Basis.
  • scalartype(d): the scalar field underlying the vector space the basis functions take value in.
  • refspace(basis): returns the ReferenceSpace of local shape functions on which the Basis is built.
  • assemblydata(basis): assemblydata returns an iterable collection elements of geometric elements and a look table ad for use in assembly of interaction matrices. In particular, for an index element_idx into elements and an index local_shape_idx in basis of local shape functions refspace(basis), ad[element_idx, local_shape_idx] returns the iterable collection of (global_idx, weight) tuples such that the local shape function at local_shape_idx defined on the element at element_idx contributes to the basis function at global_idx with a weight of weight.
  • geometry(basis): returns an iterable collection of Elements. The order in which these Elements are encountered corresponds to the indices used in the assembly data structure.

Reference Space

The reference space concept defines an API for working with spaces of local shape functions. The main role of objects implementing this concept is to allow specialization of the functions that depend on the precise reference space used.

The functions that depend on the type and value of arguments modeling reference space are:

Kernel

A kernel is a fairly simple concept that mainly exists as part of the definition of a Discrete Operator. A kernel should obey the following semantics:

In many function definitions the kernel object is referenced by operator or something similar. This is a misleading name as an operator definition should always be accompanied by the domain and range space.

Discrete Operator

Informally speaking, a Discrete Operator is a concept that allows for the computation of an interaction matrix. It is a kernel together with a test and trial basis. A Discrete Operator can be passed to assemble and friends to compute its matrix representation.

A discrete operator is a triple (kernel, test_basis, trial_basis), where kernel is a Kernel, and test_basis and trial_basis are Bases. In addition, the following expressions should be implemented and behave according to the correct semantics:

In the context of fast methods such as the Fast Multipole Method other algorithms on Discrete Operators will typically be defined to compute matrix vector products. These algorithms do not explicitly compute and store the interaction matrix (this would lead to unacceptable computational and memory complexity).

BEAST.elements — Function

elements(geo)

Create an iterable collection of the elements stored in geo. The order in which this collection produces the elements determines the index used for lookup in the data structures returned by assemblydata and quaddata.

source
BEAST.scalartype — Function
scalartype(x)

The scalar field over which the values of a global or local basis function, or an operator are defined. This should always be a scalar type, even if the basis or operator takes on values in a vector or tensor space. This data type is used to determine the eltype of assembled discrete operators.

source
BEAST.assemblydata — Function
charts, admap, act_to_global = assemblydata(basis; onlyactives=true)

Given a basis this function returns a data structure containing the information required for matrix assemble, that is, the vector charts containing Simplex elements, a variable admap of type AssemblyData, and a mapping from indices of actively used simplices to global simplices.

When onlyactives is true, another layer of indices is introduced to filter out all cells of the mesh that are not in the union of the support of the basis functions (i.e., when the basis functions are defined only on a part of the mesh).

admap is, in essence, a three-dimensional array of named tuples, which, by wrapping it in the struct AssemblyData, allows the definition of iterators. The tuple consists of the two entries

admap[i,r,c].globalindex
+admap[i,r,c].coefficient

Here, c and r are indices in the iterable set of (active) simplices and the set of shape functions on each cell/simplex: r ranges from 1 to the number of shape functions on a cell/simplex, c ranges from 1 to the number of active simplices, and i ranges from 1 to the number of maximal number of basis functions, where any of the shape functions contributes to.

For example, for continuous piecewise linear lagrange functions (c0d1), each of the three shape functions on a triangle are associated with exactly one Lagrange function, and therefore i is limited to 1.

Note: When onlyactives=false, the indices c correspond to the position of the corresponding cell/simplex whilst iterating over geometry(basis). When onlyactives=true, then act_to_global(c) correspond to the position of the corresponding cell/simplex whilst iterating over geometry(basis).

For a triplet (i,r,c), globalindex is the index in the basis of the ith basis function that has a contribution from shape function r on (active) cell/simplex c. coefficient is the coefficient of that contribution in the linear combination defining that basis function in terms of shape function.

source
BEAST.geometry — Function
geometry(basis)

Returns an iterable collection of geometric elements on which the functions in basis are defined. The order the elements are encountered needs correspond to the element indices used in the data structure returned by assemblydata.

source
BEAST.refspace — Function
refspace(basis)

Returns the ReferenceSpace of local shape functions on which the basis is built.

source
BEAST.quaddata — Function
quaddata(operator, test_refspace, trial_refspace, test_elements, trial_elements)

Returns an object cashing data required for the computation of boundary element interactions. It is up to the client programmer to decide what (if any) data is cached. For double numberical quadrature, storing the integration points for example can significantly speed up matrix assembly.

  • operator is an integration kernel.
  • test_refspace and trial_refspace are reference space objects. quadata

is typically overloaded on the type of these local spaces of shape functions. (See the implementation in maxwell.jl for an example).

  • test_elements and trial_elements are iterable collections of the geometric

elements on which the finite element space are defined. These are provided to allow computation of the actual integrations points - as opposed to only their coordinates.

source
BEAST.quadrule — Function
quadrule(operator, test_refspace, trial_refspace, test_index, test_chart, trial_index, trial_chart, quad_data)

Based on the operator kernel and the test and trial elements, this function builds an object whose type and data fields specify the quadrature rule that needs to be used to accurately compute the interaction integrals. The quad_data object created by quaddata is passed to allow reuse of any precomputed data such as quadrature points and weights, geometric quantities, etc.

The type of the returned quadrature rule will help in deciding which method of momintegrals to dispatch to.

source
BEAST.momintegrals! — Function

momintegrals!(biop, tshs, bshs, tcell, bcell, interactions, strat)

Function for the computation of moment integrals using simple double quadrature.

source
diff --git a/dev/quadstrat/index.html b/dev/quadstrat/index.html index 6dc5c97f..4da9279c 100644 --- a/dev/quadstrat/index.html +++ b/dev/quadstrat/index.html @@ -48,4 +48,4 @@ function momintegrals(op, tel, bel, qr::HighPrecisionQR) ... -end +end diff --git a/dev/tdefie/index.html b/dev/tdefie/index.html index 64a7689c..0257ae0a 100644 --- a/dev/tdefie/index.html +++ b/dev/tdefie/index.html @@ -33,4 +33,4 @@ mat"clf" patch(geo, real.(norm.(fcr))) mat"cd($(pwd()))" -mat"print('current.png', '-dpng')" +mat"print('current.png', '-dpng')" diff --git a/dev/tutorial/index.html b/dev/tutorial/index.html index db35d409..418680ef 100644 --- a/dev/tutorial/index.html +++ b/dev/tutorial/index.html @@ -16,4 +16,4 @@ efie = @discretise t[k,j]==e[k] j∈X k∈X

Solving and computing the values of the induced current in the centers of the triangles of the mesh is now straightforward:

u = solve(efie)
 fcr, geo = facecurrents(u,X)
dots out of 10: ..........
 Converting system to Matrix...done.
-LU solution of the linear system...done.

The resulting current distribution can be visualised by e.g. Matlab, Paraview, Plotly,...

+LU solution of the linear system...done.

The resulting current distribution can be visualised by e.g. Matlab, Paraview, Plotly,...