Skip to content
Snippets Groups Projects
Commit cc5c85ce authored by Alexander Müller's avatar Alexander Müller Committed by AlexanderMueller
Browse files

refactor blockvector bindings to allow for generic types

parent 7e400e26
Branches
Tags
1 merge request!562refactor blockvector bindings to allow for generic types
Pipeline #70182 passed
......@@ -111,8 +111,6 @@ namespace Dune
typedef typename BlockVector::block_type block_type;
typedef typename BlockVector::size_type size_type;
registerFieldVecMat<block_type>::apply();
using pybind11::operator""_a;
cls.def( "assign", [] ( BlockVector &self, const BlockVector &x ) { self = x; }, "x"_a );
......
......@@ -99,17 +99,28 @@ for i in range(5):
y1, y2 = blockVector(5), blockVector(5)
mat.mv(x, y1)
mat.asLinearOperator().apply(x, y2)
if (y1 - y2).two_norm > 1e-12:
raise Exception("mat.mv != mat.asLinearOperator().apply")
S = CGSolver(mat.asLinearOperator(), SeqJacobi(mat), 1e-10)
mat = None
z = blockVector(5)
_, _, converged, _, _ = S(z, y1)
if not converged:
raise Exception("CGSolver has not converged")
if (z - x).two_norm > 1e-8:
raise Exception("CGSolver unable to solve identity")
vec = dune.common.FieldVector([1])
z2 = BlockVector(vec,5)
S2=CGSolver(mat.asLinearOperator(),SeqJacobi(mat),1e-10)
_, _, converged2, _, _ = S2(z2, y2)
if not converged2:
raise Exception("CGSolver has not converged")
if (z2 - x).two_norm > 1e-8:
raise Exception("CGSolver unable to solve identity")
s = "(" + ", ".join(str(v) for v in x) + ")"
......
......@@ -55,19 +55,54 @@ def bcrsMatrix(size=0, *args, **kwargs):
return BCRSMatrix(blockType)(size,*args)
return BCRSMatrix(blockType)()
def BlockVector(typeOrBlockSize,size=None):
"""
Create a BlockVector object with specified type or block size for a FieldVector.
def BlockVector(blockSize):
if blockSize == 1:
return BlockVector1
typeName = "Dune::BlockVector< Dune::FieldVector< double ," + str(blockSize) + " > >"
Args:
typeOrBlockSize (int or str or object): The type of the block vector or the block size.
- If int, specifies the block size directly.
- If typeOrBlockSize is not an int, it should have a 'cppTypeName' attribute specifying the type and 'cppIncludes' attribute specifying additional include files.
size (int, optional): The size of the vector. Default is None.
Returns:
BlockVector: The created BlockVector object.
"""
typeName=[]
includes=["dune/python/istl/bvector.hh"]
# todo: provide other constructors
return loadvec(includes, typeName).BlockVector
if hasattr (typeOrBlockSize, "cppTypeName") and hasattr (typeOrBlockSize, "cppIncludes"):
def blockVector(size, blockSize=1):
if blockSize == 1:
typeName = f"Dune::BlockVector< {typeOrBlockSize.cppTypeName} >"
includes += typeOrBlockSize.cppIncludes
if size is None:
return loadvec(includes,typeName).BlockVector
bvec = loadvec(includes,typeName).BlockVector(size)
for i in range(len(bvec)):
bvec[i]=typeOrBlockSize
return bvec
else:
if typeOrBlockSize == 1:
if size is None:
return BlockVector1
else:
return BlockVector1(size)
typeName = "Dune::BlockVector< Dune::FieldVector< double ," + str(blockSize) + " > >"
includes = ["dune/istl/bvector.hh"]
# todo: provide other constructors
dune.common.FieldVector(typeOrBlockSize) # make sure that the FieldVector is registered
typeName = "Dune::BlockVector<Dune::FieldVector< double, " + str(typeOrBlockSize) + " > >"
if size is None:
return loadvec(includes, typeName).BlockVector
else:
return loadvec(includes, typeName).BlockVector(size)
def blockVector(size, blockSize=1):
"""
Creates a Dune BlockVector object of size `size` with FieldVector<double,`blockSize`> as blocks.
Args:
size (int): The size of the vector.
blockSize (int, optional): The block size. Default is 1.
Returns:
BlockVector: The created BlockVector object.
"""
return BlockVector(blockSize,size)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment