dune-geometry merge requestshttps://gitlab.dune-project.org/core/dune-geometry/-/merge_requests2022-05-29T16:51:42Zhttps://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/192WIP: Add a Mapping of a Geometry2022-05-29T16:51:42ZSimon PraetoriusWIP: Add a Mapping of a Geometry## Summary
This MR adds two transformations/mappings of a geometry:
1. a geometry mapped by a function
2. a geometry mapped by another geometry (chaining of two geometry)
## Examples
In the following three examples, functions defined on...## Summary
This MR adds two transformations/mappings of a geometry:
1. a geometry mapped by a function
2. a geometry mapped by another geometry (chaining of two geometry)
## Examples
In the following three examples, functions defined on different domains with different interpretations of derivatives are shown.
1. In **Example 1**, a function that can be evaluated in global coordinates with a derivative w.r.t. global coordinates is used. This requires that we map first from the reference element to the grid element and also requires that the jacobian is transformed from global coordinates to local coordinates
2. In **Example 2**, the function is defined purely local, in terms of a local basis and a linear-combination of the basis functions. Thus, the evaluation is in local coordinates and the derivatives are already w.r.t. local coordinates. No additional transformation into a real element is required and we can use the reference-element geometry instead.
3. The last example, **Example 3**, is in between the first two examples. Thereby, the function can be evaluated in local coordinates, but the derivatives are defined w.r.t. global coordinates. Thus, we need to transform only the derivative, using a special wrapper around the reference-element geometry, returning the `jacobianTransposed()` from the real geometry.
### Example 1
An analytic function for the parametrization:
```cpp
// construct a base grid e.g. with piecewise flat elements
YaspGrid<2> grid({1.0,1.0}, {10,10});
using Domain = FieldVector<double,2>;
using Range = FieldVector<double,3>;
auto sig = Functions::SignatureTag<Range(Domain)>{};
auto f = Functions::makeDifferentiableFunctionFromCallables(sig,
[](const Domain& x) -> FieldVector<double,3> {
return {x[0], x[1], std::sin(x[0] * x[1])};
},
[](const Domain& x) -> FieldMatrix<double,3,2> {
return {
{1.0, 0.0},
{0.0, 1.0},
{x[1] * std::cos(x[0] * x[1]), x[0] * std::cos(x[0] * x[1])}
};
});
for (const auto& e : elements(grid.leafGridView()))
{
// construct the mapped geometry
MappedGeometry geometry{f, e.geometry()};
}
```
#### Examples 2
A discrete function build from dune-localfunction local finite-element basis functions:
```cpp
template <class LocalBasis, class Coefficients>
class LocalFunction
{
// implementation of a discrete function based on linear-combination of basis functions
};
using LocalFiniteElement = LagrangeCubeLocalFiniteElement<double,double,2,order>;
using Range = FieldVector<double,3>;
// construct a base grid e.g. with piecewise flat elements
YaspGrid<2> grid({1.0,1.0}, {10,10});
LocalFiniteElement localFE{};
for (const auto& e : elements(grid.leafGridView()))
{
auto X = [geo=e.geometry()](const auto& local) -> Range
{
auto x = geo.global(local);
return {x[0], x[1], std::sin(x[0] * x[1])};
};
std::vector<Range> vertices;
localFE.localInterpolation().interpolate(X, vertices);
// construct the mapped geometry
LocalFunction lf{localFE.localBasis(), vertices};
MappedGeometry geometry{lf, ReferenceElementGeometry{referenceElement(e)}};
// or
// auto refElem = referenceElement(e);
// MappendGeometry geometry{lf, refElem.template geometry<0>(0)};
}
```
### Example 3
In case you get a local function created from a grid-function with a derivative in global coordinates, but the argument is already transformed by a geometry, you can use the additional wrapper `LocalDerivativeGeometry`:
```cpp
// construct a base grid e.g. with piecewise flat elements
YaspGrid<2> grid({1.0,1.0}, {10,10});
using Domain = FieldVector<double,2>;
using Range = FieldVector<double,3>;
auto sig = Functions::SignatureTag<Range(Domain)>{};
auto f = Functions::makeDifferentiableFunctionFromCallables(sig,
[](const Domain& x) -> FieldVector<double,3> {
return {x[0], x[1], std::sin(x[0] * x[1])};
},
[](const Domain& x) -> FieldMatrix<double,3,2> {
return {
{1.0, 0.0},
{0.0, 1.0},
{x[1] * std::cos(x[0] * x[1]), x[0] * std::cos(x[0] * x[1])}
};
});
// create a differentiable analytic grid function
auto fGridFunction = Functions::makeAnalyticGridViewFunction(f, grid.leafGridView());
// build the corresponding local function
auto fLocalFunction = localFunction(fGridFunction);
for (const auto& e : elements(grid.leafGridView()))
{
// bind the local function to the grid element
fLocalFunction.bind(e);
// construct the mapped geometry
MappedGeometry geometry{fLocalFunction, LocalDerivativeGeometry{e.geometry()}};
// unbind the local function from the element
fLocalFunction.unbind();
}
```https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/191Remove transition code for value semantics of ReferenceElement2022-05-29T20:40:02ZChristoph GrüningerRemove transition code for value semantics of ReferenceElementDeprecation message for Transitional::ReferenceElement not working for Clang (due to a [Clang bug](https://github.com/llvm/llvm-project/issues/18236)) and not working for dune-grid's GeometryGrid (don't know why).Deprecation message for Transitional::ReferenceElement not working for Clang (due to a [Clang bug](https://github.com/llvm/llvm-project/issues/18236)) and not working for dune-grid's GeometryGrid (don't know why).https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/190WIP: Create DUNE_PROJECT_TARGET2022-02-09T23:02:13ZSimon PraetoriusWIP: Create DUNE_PROJECT_TARGETThis MR is associated to the core/dune-common!944This MR is associated to the core/dune-common!944https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/189WIP: [python] Without Python bindings, show Python tests as skipped2022-01-26T23:19:09ZChristoph GrüningerWIP: [python] Without Python bindings, show Python tests as skippedhttps://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/175Wrap the quadrature floating-point numbers in the macro DUNE_NUMBER2021-11-20T12:05:57ZSimon PraetoriusWrap the quadrature floating-point numbers in the macro DUNE_NUMBER### Summary
The maxima script to generate 1d quadrature rules produces two copies of the same code: one with floating point numbers, one with strings. While it is generated automatically, it is a lot of code duplication that is actually ...### Summary
The maxima script to generate 1d quadrature rules produces two copies of the same code: one with floating point numbers, one with strings. While it is generated automatically, it is a lot of code duplication that is actually not necessary. Also for future high-precision rules, duplicating all the code is not a proper solution. This MR introduces a small helper macro that expands the argument to both, a floating point constant and a string. Additionally a helper function is introduced to convert from the two generated arguments to the output type.
### Discussion
The floating-point number type is currently fixed to `double`. Actually, one could pass the numbers as `long double` by adding the literal suffix `L`. Unfortunately, the `GMPField` type does not allow to be constructed from string representations of numbers that include any suffix. Thus, one would have to workaround this, e.g., by stripping off any suffix before passing to the constructor of the number type.https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/174Use a container with flexible key type as cache for quadrature rules2021-07-04T14:26:06ZSimon PraetoriusUse a container with flexible key type as cache for quadrature rules### Summary
For a thread-safe cache of quadrature rules, currently a vector of vectors of vectors with several `std::call_once` flags is used. This makes reading the code very complicated and does not necessarily lead to any better perfo...### Summary
For a thread-safe cache of quadrature rules, currently a vector of vectors of vectors with several `std::call_once` flags is used. This makes reading the code very complicated and does not necessarily lead to any better performance than the proposed solution in this MR. Replace the nested vector with a simple `std::map` where the key allows to identify all quadrature rules and additionally makes it possible to extend quadrature rules properties in the future.
### Discussion
1. Instead of a `std::map` container, any other associative container could be used. I have chose a simple map since comparison of three integers by `<` is very simple to implement. The alternative `std::unordered_map` would require a hash function. While this is possible, it requires more code. I do not expect a measurable performance difference.
2. The cache is implemented as `thread_local`. This allows to omit more complicated locking mechanisms and should be thread-safe by definition. It is, however, unclear to me, whether this has negative performance consequences it threads are closed and reopened multiple times, i.e., if this makes it necessary to re-"compute"/re-fill the cache or whether this preserves the previously cached values. An alternative implementation uses the `static` keyword:
```c++
// Container to store the cached values
static std::map<QuadratureKey, QuadratureRule> quadCache;
// mutex used to access the data in the container, necessary since
// access emplace is read-write.
using mutex_type = std::shared_timed_mutex;
static mutex_type access_mutex;
// define the quadrature key to be used to access the quadrature rule
const QuadratureKey key{t.id(),p,qt};
// first try to lock for read-only, if a quadratur rule for key is found, return it,
// if not, obtain a unique_lock to insert a new rule.
std::shared_lock<mutex_type> read_lock(access_mutex);
auto it = quadCache.find(key);
if (it != quadCache.end())
return it->second;
else {
read_lock.unlock();
QuadratureRule rule = QuadratureRuleFactory<ctype,dim>::rule(t,p,qt);
std::unique_lock<mutex_type> write_lock(access_mutex);
auto new_it = quadCache.emplace(key, std::move(rule));
return new_it.first->second;
}
```https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/173Make the enum QuadratureTypes::Enum a proper enum class2021-08-08T09:17:44ZSimon PraetoriusMake the enum QuadratureTypes::Enum a proper enum class### Summary
Since c++11 we have enum classes to solve the problem of scoped names for enums. This MR replaces the old workaround with a proper enum class for `QuadratureType`. This will allow in future c++ standards to `use QuadratureTyp...### Summary
Since c++11 we have enum classes to solve the problem of scoped names for enums. This MR replaces the old workaround with a proper enum class for `QuadratureType`. This will allow in future c++ standards to `use QuadratureTypes` inside of switch case statements and makes the code more readable.
### Discussion
In order to use enum class types as integers, e.g. as index in a vector, it must be explicitly converted into the underlying integer type. No implicit conversion happens. Once the quadrature cache is replaced with a proper map, this cast is not longer necessary. See !174.https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/170Add geometry mapped by an element-function2022-02-25T17:49:08ZSimon PraetoriusAdd geometry mapped by an element-function## Summary
A geometry that is defined as the mapping of a reference geometry
#### Details
A geometry can be seen as a class interface representing a reference domain `R` plus a mapping `f`, such that `f(R)` is of the same topology as `R...## Summary
A geometry that is defined as the mapping of a reference geometry
#### Details
A geometry can be seen as a class interface representing a reference domain `R` plus a mapping `f`, such that `f(R)` is of the same topology as `R` and `f` a bijective (diffeomorphic) mapping between `R` and `f(R)`. The `MappedGeometry` provides this interface for a given `ReferenceElement` `R` and mapping `f`.
Additionally, the `MappedGeometry` allows to be defined on a different reference domain than associated to the mapping `f`. This could be, for example, the domain of a subEntity or element-intersection. In this case, we build a composition of two mappings `f` and `g`, where `g` represents a local-geometry transform between the two reference domains. Let `I` be the reference domain of the subEntity and `g` the associated geometry mapping (that is bijective between `I` and `g(I)`), then the `MappedGeometry` represents the composition `f o g` associated to the reference domain `I`.
We require for the local function `f` to fulfill some interface to be allowed in this wrapper class. Let `lf` be an object of local function, `x` be a local coordinate in the reference element of the geometry, then (formulated in c++-20 concepts)
```c++
template <class LF, class Range, class Domain, class LocalContext,
class DLF = std::decay_t<decltype(derivative(std::declval<const LF&>())> >
concept LocalFunction = require (const LF& lf, const Domain& x)
{
{ lf(x) } -> std::convertible_to<Range>;
{ lf.localContext() } -> std::convertible_to<LocalContext>;
{ derivative(lf) } -> std::regular_invocable<Domain>;
require (const DLF& dlf, const Domain& x)
{
{ dlf(x) } -> std::convertible_to<typename DerivativeTraits<Range(Domain)>::Range>;
};
};
```
In this concept we make explicit notion of a "LocalContext" the local-function is associated to. This context refers to an entity or intersection the local-function was bound. This means that the local-function can be evaluated in the local-coordinates of that `LocalContext`. We thus require, that the local-context provides also a geometry.
## Example 1
This example uses *dune-functions*, just to make the example short. One can easily define a local function following the `test-mappedgeometry.cc`
```c++
// construct a base grid e.g. with piecewise flat elements
YaspGrid<2> grid({1.0,1.0}, {10,10});
using Domain = FieldVector<double,2>;
using Range = FieldVector<double,3>;
auto sig = Functions::SignatureTag<Range(Domain)>{};
auto f = Functions::makeDifferentiableFunctionFromCallables(sig,
[](const Domain& x) -> FieldVector<double,3> {
return {x[0], x[1], std::sin(x[0] * x[1])};
},
[](const Domain& x) -> FieldMatrix<double,3,2> {
return {
{1.0, 0.0},
{0.0, 1.0},
{x[1] * std::cos(x[0] * x[1]), x[0] * std::cos(x[0] * x[1])}
};
});
// create a differentiable analytic grid function
auto fGridFunction = Functions::makeAnalyticGridViewFunction(f, grid.leafGridView());
// build the corresponding local function
auto fLocalFunction = localFunction(fGridFunction);
for (const auto& e : elements(grid.leafGridView()))
{
// bind the local function to the grid element
fLocalFunction.bind(e);
// construct the mapped geometry
MappedGeometry geometry{referenceElement(e), fLocalFunction};
// unbind the local function from the element
fLocalFunction.unbind();
}
```
## Example 2
In this example also intersection geometries are transformed. Since the local function must be bound to a grid element, we need an additional geometry transformation from the intersection geometry to the element geometry. The local-geometry transform is passed as additional argument. It can also be used to define a geometry on sub-entities of the grid element, by using a reference sub-element geometry mapping.
Let `grid`, `f`, `fGridFunction`, `fLocalFunction` be defined as in Example 1.
```c++
for (const auto& e : elements(grid.leafGridView()))
{
// bind the local function to the grid element
fLocalFunction.bind(e);
for (const auto& is : intersections(element, grid.leafGridView())
{
// construct the mapped geometry for the intersection
MappedGeometry isGeometry{referenceElement(is), fLocalFunction, is.geometryInInside()};
}
auto refElem = referenceElement(e);
for (int i = 0; i < refElem.size(1); ++i)
{
auto edgeRefElem = referenceElement<double,grid.dimension-1>(refElem.type(i,1));
// construct a mapped geometry for the i'th subentity of codim 1
// NOTE: this required a well-oriented (twist-free) grid
MappedGeometry edgeGeometry{edgeRefElem, fLocalFunction, refElem.geometry<1>(i)};
}
// unbind the local function from the element
fLocalFunction.unbind();
}
```https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/169Add parametrized (curved) geometry2022-02-25T18:17:22ZSimon PraetoriusAdd parametrized (curved) geometry### Summary
This is the first version of a parametrized geometry that allows for curved elements. The parametrization is based on a local finite-element map. Note, the local finite-element is just a template parameter and thus not a depe...### Summary
This is the first version of a parametrized geometry that allows for curved elements. The parametrization is based on a local finite-element map. Note, the local finite-element is just a template parameter and thus not a dependency to dune-localfunctions.
#### Example 1
Parametrize a geometry by 2nd order Lagrange elements:
```c++
// construct a base grid e.g. with piecewise flat elements
YaspGrid<2> grid({1.0,1.0}, {10,10});
// geometry polynomial order
const int order = 2;
// define a local finite-element (from dune-localfunction)
using LocalFiniteElement = LagrangeCubeLocalFiniteElement<double,double,2,order>;
LocalFiniteElement localFE{};
for (const auto& e : elements(grid.leafGridView()))
{
// projection from local coordinates in the reference element
// to 2d coordinates in the base grid to 3d global coordinates
auto X = [geo=e.geometry()](const auto& local) -> FieldVector<double,3>
{
auto x = geo.global(local);
return {x[0], x[1], std::sin(x[0] * x[1])};
};
// construct the parametrized geometry
ParametrizedGeometry geometry{referenceElement(e), localFE, X};
}
```
#### Example 2
First construct a global grid parametrization in form of a vector holding the coordinates of all Lagrange points, then locally build the corresponding geometry. Note, the example uses a global basis from dune-functions, but it can be implemented also low-level
```c++
// construct a base grid e.g. with piecewise flat elements
YaspGrid<2> grid({1.0,1.0}, {10,10});
using namespace Dune::Functions::BasisFactory;
auto parametrizationBasis
= makeBasis(grid.leafGridView(), power<3>(lagrange<2>(), blockedInterleaved()));
// create avector that holds all the Lagrange points on the curved grid
std::vector<FieldVector<double,3>> points(parametrizationBasis.size());
Functions::interpolate(parametrizationBasis, points,
[](auto const& x) -> FieldVector<double,3> {
return {x[0], x[1], std::sin(x[0] * x[1])};
});
auto localView = parametrizationBasis.localView();
for (const auto& e : elements(grid.leafGridView()))
{
localView.bind(e);
auto const& node = localView.tree().child(0);
// extract the points on the current elements
std::vector<FieldVector<double,3>> localPoints(node.size());
for (std::size_t i = 0; i < node.size(); ++i) {
std::size_t idx = localView.index(node.localIndex(i))[0];
localPoints[i] = points[idx];
}
// construct the parametrized geometry
ParametrizedGeometry geometry{referenceElement(e), node.finiteElement(), localPoints};
}
```https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/155WIP: fix pkg-config2021-01-19T11:24:14ZChristian EngwerWIP: fix pkg-configCorresponding to core/dune-common!910 this MR is a start to fix `pkg-config` support.Corresponding to core/dune-common!910 this MR is a start to fix `pkg-config` support.https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/153[feature][Python] Add default reference-element containers to the library.2022-01-27T12:55:56ZRobert K[feature][Python] Add default reference-element containers to the library.This is done for dim <= 3 to reduce compile times.This is done for dim <= 3 to reduce compile times.https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/144use pip installed pybind11 and compile _geometry.so in dune-py module2022-01-09T15:56:00ZAndreas Dedneruse pip installed pybind11 and compile _geometry.so in dune-py modulehttps://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/122Set upper bounds on dune module dependencies2019-03-07T13:16:48ZLinus SeelingerSet upper bounds on dune module dependenciesEnsure no incompatible future releases of dependencies are mixed with this module.Ensure no incompatible future releases of dependencies are mixed with this module.https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/105WIP: Add StaticGeometryType2019-01-30T16:48:27ZCarsten GräserWIP: Add StaticGeometryTypeThis encodes all the information of a `GeometryType` as template
parameters. Some remarks on the implementation:
* This may be a solition for using `GeometryType`s a template paremters.
* All methods in `Dune::GeometryTypes::` have analo...This encodes all the information of a `GeometryType` as template
parameters. Some remarks on the implementation:
* This may be a solition for using `GeometryType`s a template paremters.
* All methods in `Dune::GeometryTypes::` have analogouges in
`Dune::StaticGeometryTypes::` and additional template
aliases for parametrized types. Hence you can e.g. use Simplex<dim>.
* There is `LocalHybridGeometryTypeIndex` which works like `LocalGeometryTypeIndex`
but supports GeometryType and StaticGeometryType.
* Some `topologyIds` are different: If you want to overload for specific
`StaticGeometryType`s, you must be sure that the same geometry type
leads to the same type. Hence `topologyIds` are normalize to have 0
in the first bit.
* One may consider using `uniqueTopologyId = (topologyId >> 1)`
as template parameter instead to get rid of the nasty non-uniqueness.
For now this is a **non-intrusive proposal for discussion** and thus marked WIP.
One could also merge `Dune::StaticGeometryTypes` into `Dune::GeometryTypes::` and
`LocalHybridGeometryTypeIndex` into `LocalGeometryTypeIndex`.
This would change some dynamic types to static ones.
However, this should not be a problem since each `StaticGeometryType`
can be converted into its dynamic counterpart.https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/103Draft: Step 2 of redesign of quadrature rules2021-07-04T13:53:40ZSimon PraetoriusDraft: Step 2 of redesign of quadrature rules### Summary
The classes `QuadratureRule`, `QuadratureRules` and `QuadratureRuleFactory` are cleaned up in style, using static map for storing rules and modern enum classes. The quadrature rules are put into a common programming style, e....### Summary
The classes `QuadratureRule`, `QuadratureRules` and `QuadratureRuleFactory` are cleaned up in style, using static map for storing rules and modern enum classes. The quadrature rules are put into a common programming style, e.g. no intermediate point/weight construction, but directly storing coefficients in the quadrature rule.
### Details
Four major changes are introduced in this MR:
1. Quadrature rules (points and weights) are constructed directly and not via an intermediate c-array.
2. The quadrature rules are stored in a cache in form of a `thread_local std::map`. This allows to easily extend the key types (currently to topology-id, the quadrature order and quadrature type). In some initial tests there was no performance different measurable. It could be discussed whether it would be better to use a `static std::map` + locking, or whether another map, like an `std::unordered_map` would be the better choice.
3. The `QuadratureType` enum is now an `enum class`. This might be a breaking change in some use-code, if the quadrature type is stored somewhere. Needed changes in user code: replace `Dune::QuadratureType::Enum` by `Dune::QuadratureType`.
4. Instead of repeating all quarature points and weights for flating-point numbers and for string construction, a utility function/macro `DUNE_NUMBER -> Dune::Impl::number` is introduced to cover both cases at once. The macro takes a type and a floating-point literal and expands the latter into a `double` and a `const char*`. In the `Dune::Impl::number` template, it is then decided whether to construct from a string or returning the double value directly.
### Motivation
This is the **second** step of a cleanup and redesign of the quadrature rules: 1.) put all rules into a similar structure. 2.) cleanup of the classes QuadratureRules and QuadratureFactory, 3.) Automatically generate coefficients for all rules with higher precision
The implementation intends to establish a common style that can be generated automatically in step 3 from a database.
### Remarks
The simplexquadraturerule is transform from the old code using the attached python [convert_rule.py](/uploads/52f8b9d3e1f0fc90759e6c58d6dd79ef/convert_rule.py) and some postprocessing.https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/69WIP: caching axisalignedcubegeometry2017-10-05T12:42:43ZDominic Kempfdominic.kempf@iwr.uni-heidelberg.deWIP: caching axisalignedcubegeometryThe branch implements a caching version of `AxisAlignedCubeGeometry`. It pre-evaluates
and stores geometric quantities such as Jacobian(Inverse)Transposed, corners and volume
upon geometry instantiation. The implementation is switched on...The branch implements a caching version of `AxisAlignedCubeGeometry`. It pre-evaluates
and stores geometric quantities such as Jacobian(Inverse)Transposed, corners and volume
upon geometry instantiation. The implementation is switched on through a template parameter,
that defaults to `False`.
Right now, I am using this to investigate a performance bottleneck in my application.
I hope I can gather some insight on if and when this implementation is to be favored.
Not opting to merge right now, I opened the MR for visibility and feedback if someone has
been toying around with that before.https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/30[wip] Use reference qualifiers to avoid reference to temporaries2016-10-12T09:57:33ZCarsten Gräser[wip] Use reference qualifiers to avoid reference to temporariesUntil now one could write code like
auto&& jt = AffineGeometry<...>(...).jacobianTransposed(...);
This compiles and (depending on the surrounding code)
may even lead to the expected results. However, since
AffineGeometry hands out a ...Until now one could write code like
auto&& jt = AffineGeometry<...>(...).jacobianTransposed(...);
This compiles and (depending on the surrounding code)
may even lead to the expected results. However, since
AffineGeometry hands out a reference to the stored
matrix jt is a reference to a temporary that may
be overwritten by other stack variables.
By explicitly deleting the overload for r-values the compiler
rejects the above code and only allows to call jacobianTransposed()
on l-values.https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/13[WIP] Allow vectorization types to be used with dune-geometry2020-02-22T11:15:48ZJö Fahlkejorrit@jorrit.de[WIP] Allow vectorization types to be used with dune-geometryThis branch builds on core/dune-common!81. It makes it possible to stuff vectorized types from some vectorization library into the templates provided by dune-geometry.
See also: core/dune-common!81, core/dune-istl!17.
[WIP] Check ...This branch builds on core/dune-common!81. It makes it possible to stuff vectorized types from some vectorization library into the templates provided by dune-geometry.
See also: core/dune-common!81, core/dune-istl!17.
[WIP] Check that it still works with current dune-common (i.e. do a unit test).DUNE 3.0.0https://gitlab.dune-project.org/core/dune-geometry/-/merge_requests/4Feature/#778 failing pseudo inverse2017-07-31T17:36:29ZElias PippingFeature/#778 failing pseudo inverseThis is an attempt to fix flyspray/FS#778 and flyspray/FS#866.
3x3 matrices can now be inverted via `rightInvA()` that would previously yield an error (this is not tested anywhere).
2x3 matrices can now be pseudo-inverted via `righ...This is an attempt to fix flyspray/FS#778 and flyspray/FS#866.
3x3 matrices can now be inverted via `rightInvA()` that would previously yield an error (this is not tested anywhere).
2x3 matrices can now be pseudo-inverted via `rightInvA()` that would previously yield an error (this is tested in the testmatrix test). What is the quality of said 3x2 inverse (called `invB`)? Here's the output that the test produces for me:
```
2: Running tests for type: double
2: sqrtDetAAT = 2.22045e-16
2: detA = 2.22045e-16
2: invA = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: invA*A*invA = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: Identity invA*A*invA = invA satisfied up to error: 0
2: A = [ 0.1 -0.01, 0.1 -0.01 ]
2: A*invA*A = [ 0.1 -0.01, 0.1 -0.01 ]
2: Identity A*invA*A = A satisfied up to error: 0
2: sqrtDetBBT = 2.22045e-16
2: detB = 2.22045e-16
2: invB = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: invB*B*invB = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: Identity invB*B*invB = invB satisfied up to error: 0
2: B = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: B*invB*B = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: Identity B*invB*B = B satisfied up to error: 0
2:
2: Running tests for type: long double
2: sqrtDetAAT = 2.22045e-16
2: detA = 2.22045e-16
2: invA = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: invA*A*invA = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: Identity invA*A*invA = invA satisfied up to error: 0
2: A = [ 0.1 -0.01, 0.1 -0.01 ]
2: A*invA*A = [ 0.1 -0.01, 0.1 -0.01 ]
2: Identity A*invA*A = A satisfied up to error: 0
2: sqrtDetBBT = 2.22045e-16
2: detB = 2.22045e-16
2: invB = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: invB*B*invB = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: Identity invB*B*invB = invB satisfied up to error: 0
2: B = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: B*invB*B = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: Identity B*invB*B = B satisfied up to error: 0
2:
2: Running tests for type: GMPField<72>
2: sqrtDetAAT = 2.22045e-16
2: detA = 2.22045e-16
2: invA = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: invA*A*invA = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: Identity invA*A*invA = invA satisfied up to error: 1.50463e-38
2: A = [ 0.1 -0.01, 0.1 -0.01 ]
2: A*invA*A = [ 0.1 -0.01, 0.1 -0.01 ]
2: Identity A*invA*A = A satisfied up to error: 1.61633e-41
2: sqrtDetBBT = 2.22045e-16
2: detB = 2.22045e-16
2: invB = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: invB*B*invB = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: Identity invB*B*invB = invB satisfied up to error: 1.17549e-38
2: B = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: B*invB*B = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: Identity B*invB*B = B satisfied up to error: 1.61633e-41
2:
2: Running tests for type: GMPField<160>
2: sqrtDetAAT = 2.22045e-16
2: detA = 2.22045e-16
2: invA = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: invA*A*invA = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: Identity invA*A*invA = invA satisfied up to error: 8.44269e-58
2: A = [ 0.1 -0.01, 0.1 -0.01 ]
2: A*invA*A = [ 0.1 -0.01, 0.1 -0.01 ]
2: Identity A*invA*A = A satisfied up to error: 1.20888e-60
2: sqrtDetBBT = 2.22045e-16
2: detB = 2.22045e-16
2: invB = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: invB*B*invB = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: Identity invB*B*invB = invB satisfied up to error: 6.37237e-58
2: B = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: B*invB*B = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: Identity B*invB*B = B satisfied up to error: 1.20888e-60
```
(Update: redid the 2x3 case. Relying on maple was not a good idea. The generated code was unstable and incomprehensible. Now it behaves just as intended.)
In summary: the 2x3 pseudo-inverse is perfectly usable for the test case, also with the `double` field type.Oliver SanderOliver Sander