dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #20857
Some patches for dolfin
I've attached some patches for dolfin. Please apply (or not). Short
description follows, and a question...
1) Allow tensor Constants in python (initialised from e.g. a numpy matrix)
2) Define a split method on FunctionSpace in python. This seems natural
since Function has one... but not important.
3) Define the methods "near" and "between" in both compiled and python
Expressions (etc). So you can write
def boundary(self):
return near(x[0], 0)
or
boundary = compile_subdomains("near(x[0], 0)")
I think "near" is good and gets rid of many of those shouting
DOLFIN_EPSs. "between" is perhaps more questionable... I leave it to you.
4) This matches up with a viper patch sent to fenics-viper@lp, and
simplifies plot() a bit. It also does less projection of those types viper
can actually handle --- which means that for example DG plots are slightly
different (but maybe not better?).
5) I find it easier to understand DIVERGED_INDEFINITE_MAT than -8 as an
error message.
6) PETSc development version requires some changes
7) This one is really unimportant and you may not agree anyway.
8) This one IS important. Allows assembly into rectangular matrices with
epetra. This used to crash, and tracing it down was a bitch due to epetra
throwing bare-int exceptions (which are rethrown by swig).
9) A utility function on DirichletBC for symmetric modification of block
matrices. Also implements EpetraMatrix::set_row. These can probably be made
much more efficient by someone who knows the la backends better (I threw in
a "this is inefficient" log msg there to make this point), but it seems to
work in superficial testing.
Now the question: I seem to hit dolfin asserts a lot when I try things, and
that's annoying since they call abort() and thus kill (i)python. I see that
they used to be implemented as runtime exceptions, which plays much better
with python, but that was changed. Was there a good reason for that change,
or can we bring back exceptions? (the change itself is trivial and I can
send a patch, but there was probably a reason...)
--
-j.
From c64191b65c81593675e090d789603a8ce2d10f2a Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Mon, 10 Jan 2011 10:09:24 +0100
Subject: [PATCH 1/9] Implement tensor constants in python layer.
---
site-packages/dolfin/function/constant.py | 34 +++++++++++++---------------
1 files changed, 16 insertions(+), 18 deletions(-)
diff --git a/site-packages/dolfin/function/constant.py b/site-packages/dolfin/function/constant.py
index 7886e85..9d7e169 100644
--- a/site-packages/dolfin/function/constant.py
+++ b/site-packages/dolfin/function/constant.py
@@ -11,33 +11,31 @@ __all__ = ["Constant"]
# Import UFL and SWIG-generated extension module (DOLFIN C++)
import ufl
import dolfin.cpp as cpp
+import numpy
class Constant(ufl.Coefficient, cpp.Constant):
def __init__(self, value):
"""Create constant-valued function with given value.
- The value may be either a single scalar value or a tuple/list of
- values for vector-valued functions."""
+ The value may be either a single scalar value, or a tuple/list of
+ values for vector-valued functions, or nested lists or a numpy array for
+ tensor-valued functions."""
- # Create UFL element
- if isinstance(value, (float, int)):
+ array = numpy.array(value)
+ dim = len(array.shape)
+ floats = map(float, array.flat)
+
+ # Create UFL element and initialize constant
+ if dim == 0:
self._ufl_element = ufl.FiniteElement("Discontinuous Lagrange", None, 0)
- elif isinstance(value, (tuple, list)):
- self._ufl_element = ufl.VectorElement("Discontinuous Lagrange", None, 0, len(value))
+ cpp.Constant.__init__(self, floats[0])
+ elif dim == 1:
+ self._ufl_element = ufl.VectorElement("Discontinuous Lagrange", None, 0, len(floats))
+ cpp.Constant.__init__(self, floats)
else:
- cpp.error("Unable to create Constant, unhandled type: %s" % str(value))
+ self._ufl_element = ufl.TensorElement("Discontinuous Lagrange", None, 0, shape=array.shape)
+ cpp.Constant.__init__(self, list(array.shape), floats)
# Initialize base classes
ufl.Coefficient.__init__(self, self._ufl_element)
- cpp.Constant.__init__(self, _to_float(value))
-
-def _to_float(value):
- "Convert to float or a (nested) list of floats."
-
- if isinstance(value, (float, int)):
- return float(value)
- elif isinstance(value, (tuple, list)):
- return [_to_float(v) for v in value]
- else:
- cpp.error("Unable to interpret value for constant: " + str(value))
--
1.7.1
From cf4f3b4abf365b2845dddaf1371d0d43343f982a Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Mon, 10 Jan 2011 14:35:52 +0100
Subject: [PATCH 2/9] Define a split() method on FunctionSpaces (python).
---
site-packages/dolfin/function/functionspace.py | 8 +++++++-
1 files changed, 7 insertions(+), 1 deletions(-)
diff --git a/site-packages/dolfin/function/functionspace.py b/site-packages/dolfin/function/functionspace.py
index 7533c44..d3eea68 100644
--- a/site-packages/dolfin/function/functionspace.py
+++ b/site-packages/dolfin/function/functionspace.py
@@ -97,11 +97,17 @@ class FunctionSpaceBase(cpp.FunctionSpace):
if self.num_sub_spaces() == 1:
raise ValueError, "no SubSpaces to extract"
if i >= self.num_sub_spaces():
- raise ValueError, "Can only extract SubSpaces with i = 0..%d"%self.num_sub_spaces()
+ raise ValueError, "Can only extract SubSpaces with i = 0..%d"%(self.num_sub_spaces()-1)
assert(hasattr(self._ufl_element,"sub_elements"))
element = self._ufl_element.sub_elements()[i]
return FunctionSpaceFromCpp(cpp.SubSpace(self,i), element)
+ def split(self):
+ S = []
+ for i in range(self.num_sub_spaces()):
+ S.append(self.sub(i))
+ return S
+
class FunctionSpaceFromCpp(FunctionSpaceBase):
" FunctionSpace represents a finite element function space."
def __init__(self, cppV, element=None):
--
1.7.1
From 824e413bc6842a46a44ecc9262b81cf06cd389e9 Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Wed, 12 Jan 2011 16:01:04 +0100
Subject: [PATCH 3/9] Introduce the functions "near" and "between" in header for compiled modules as well as python.
These are commonly used functions in subdomain definitions, as
"near(x[0], coord)" -> "abs(x[0]-coord) < eps" or
"between(coord1, x[0], coord2)" -> "coord0-eps < x[0] < coord1+eps".
---
site-packages/dolfin/__init__.py | 2 +-
site-packages/dolfin/compilemodules/__init__.py | 43 ++++++++++++++++++----
2 files changed, 36 insertions(+), 9 deletions(-)
diff --git a/site-packages/dolfin/__init__.py b/site-packages/dolfin/__init__.py
index d0cfc2f..b8140a9 100644
--- a/site-packages/dolfin/__init__.py
+++ b/site-packages/dolfin/__init__.py
@@ -26,7 +26,7 @@ from dolfin.common.plot import *
from dolfin.common.time import *
from dolfin.common.memory import *
-from dolfin.compilemodules import compile_extension_module
+from dolfin.compilemodules import compile_extension_module, near, between
from dolfin.compilemodules.expressions import compile_expressions
from dolfin.compilemodules.subdomains import compile_subdomains
from dolfin.compilemodules.jit import jit
diff --git a/site-packages/dolfin/compilemodules/__init__.py b/site-packages/dolfin/compilemodules/__init__.py
index b7e25a9..04b723f 100644
--- a/site-packages/dolfin/compilemodules/__init__.py
+++ b/site-packages/dolfin/compilemodules/__init__.py
@@ -20,7 +20,8 @@ from dolfin.compilemodules.jit import mpi_jit_decorator
__all__ = ["compile_extension_module",
"expression_to_code_fragments",
- "math_header"]
+ "math_header",
+ "near", "between"]
# A list of supported math builtins
_math_builtins = [
@@ -30,9 +31,26 @@ _math_builtins = [
"exp", "frexp", "ldexp", "log", "log10", "modf",
"pow", "sqrt", "ceil", "fabs", "floor", "fmod"]
-math_header = "\n// cmath functions\n" + \
- "\n".join("using std::%s;"%mf for mf in _math_builtins) + \
- "\nconst double pi = DOLFIN_PI;\n\n"
+_exclude_as_variables = ["pi", "near", "between", "true", "false"]
+
+math_header = """
+// cmath functions
+
+%s
+
+const double pi = DOLFIN_PI;
+
+static inline bool between(double x0, double x, double x1)
+{
+ return (x0-DOLFIN_EPS < x && x < x1+DOLFIN_EPS);
+}
+
+static inline bool near(double x, double x0)
+{
+ return between(x0,x,x0);
+}
+
+""" % ("\n".join("using std::%s;"%mf for mf in _math_builtins))
_cpp_keywords = ["auto","const","double","float","int","short","struct","unsigned",
"break","continue","else","for","long","signed","switch","void",
@@ -40,12 +58,21 @@ _cpp_keywords = ["auto","const","double","float","int","short","struct","unsigne
"char","do","extern","if","return","static","union","while",
"asm","dynamic_cast","namespace","reinterpret_cast","try",
"bool","explicit","new","static_cast","typeid","volatile",
- "catch","false","operator","template","typename",
+ "catch","operator","template","typename",
"class","friend","private","this","using",
"const_cast","inline","public","throw","virtual",
- "delete","mutable","protected","true","wchar_t",
+ "delete","mutable","protected","wchar_t",
"or","and","xor","not"]
+# Define these to match the C++ header functions
+epsilon = dolfin.DOLFIN_EPS
+def between(x0, x, x1):
+ """Returns true if x0-epsilon < x < x1+epsilon"""
+ return x0-epsilon < x < x1+epsilon
+def near(x0, x):
+ """Returns true if x0-epsilon < x < x0+epsilon"""
+ return x0-epsilon < x < x0+epsilon
+
def expression_to_code_fragments(expr, defaults, arguments):
"A help function which extract a dict with code snippets from an expression"
@@ -71,8 +98,8 @@ def expression_to_code_fragments(expr, defaults, arguments):
# Remove any variables defined in the arguments list
variables.difference_update(arguments)
- # Remove the builtin math functions from the variables (and 'pi', we use DOLFIN_PI)
- variables.difference_update(_math_builtins+["pi"])
+ # Remove the builtin math functions from the variables (and other definitions)
+ variables.difference_update(_math_builtins+_exclude_as_variables)
# Remove the numerals from the variables
numerals = [v for v in variables if v[0] in "0123456789"]
--
1.7.1
From 3ac4d409631decdcd45c558b9b3a6301a2a98672 Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Wed, 12 Jan 2011 17:18:38 +0100
Subject: [PATCH 4/9] Simplify plot.py, and let viper decide which types are supported.
Requires update of fenics_viper.
---
site-packages/dolfin/common/plot.py | 66 ++++++++++++++--------------------
1 files changed, 27 insertions(+), 39 deletions(-)
diff --git a/site-packages/dolfin/common/plot.py b/site-packages/dolfin/common/plot.py
index f760854..dc81002 100644
--- a/site-packages/dolfin/common/plot.py
+++ b/site-packages/dolfin/common/plot.py
@@ -17,58 +17,45 @@ import ufl
__all__ = ['Viper', 'plot', 'update', 'interactive', 'save_plot', 'figure']
-# Intelligent plot command that handles projections and different objects
-def dolfin_plot(object, *args, **kwargs):
- "Plot given object (Mesh, MeshFunction or Function)"
-
- from viper.viper_dolfin import plot as viper_plot
- from ffc import plot as ffc_plot
-
- # Plot mesh
- if isinstance(object, cpp.Mesh):
- return viper_plot(object, *args, **kwargs)
-
- # Plot mesh function
- if isinstance(object, cpp.MeshFunctionUInt) or \
- isinstance(object, cpp.MeshFunctionInt) or \
- isinstance(object, cpp.MeshFunctionDouble):
- return viper_plot(object, *args, **kwargs)
+def make_viper_object(object):
+ for plottype,classes in viper_dolfin.plottable:
+ if isinstance(object, classes):
+ return object
- # Plot function
- if isinstance(object, cpp.Function):
- return viper_plot(object, *args, **kwargs)
-
- # Plot expression
- if isinstance(object, cpp.Expression):
- if "mesh" not in kwargs:
- raise TypeError, "expected a mesh when plotting an expression."
- return viper_plot(object, *args, **kwargs)
-
- # Plot function plot data
- if isinstance(object, cpp.FunctionPlotData):
- return viper_plot(object, *args, **kwargs)
-
- # Plot boundary condition
if isinstance(object, cpp.DirichletBC):
bc = object
V = bc.function_space()
v = Function(V)
bc.apply(v.vector())
- return viper_plot(v, *args, **kwargs)
-
- # Plot element
- if isinstance(object, ufl.FiniteElementBase):
- return ffc_plot(object, *args, **kwargs)
+ return v
# Try projecting function
from dolfin.fem.project import project
try:
- print "Object cannot be plotted directly, projecting to piecewise linears."
u = project(object)
- return plot(u, *args, **kwargs)
+ cpp.info("Object cannot be plotted directly, projecting to piecewise linears.")
+ return u
except:
raise RuntimeError, ("Don't know how to plot given object and projection failed: " + str(object))
+
+# Intelligent plot command that handles projections and different objects
+# (aliased as 'plot')
+def dolfin_plot(object, *args, **kwargs):
+ "Plot given object (Mesh, MeshFunction or Function)"
+
+ # Plot element
+ if isinstance(object, ufl.FiniteElementBase):
+ import ffc
+ return ffc.plot(object, *args, **kwargs)
+
+ # Check expression
+ if isinstance(object, cpp.Expression) and "mesh" not in kwargs:
+ raise TypeError, "expected a mesh when plotting an expression."
+
+ return viper_dolfin.plot(make_viper_object(object), *args, **kwargs)
+
+
# Check DOLFIN_NOPLOT
do_nothing = False
if os.environ.has_key('DOLFIN_NOPLOT'):
@@ -78,11 +65,12 @@ else:
# Check for Viper
try:
+ from viper import viper_dolfin
for x in __all__:
exec ('from viper.viper_dolfin import %s' % x)
plot = dolfin_plot
except ImportError, details:
- cpp.info(str(details))
+ cpp.warning(str(details))
cpp.warning("Unable to import Viper, plotting disabled.")
do_nothing= True
--
1.7.1
From c7e68ebaf7a366ab930200a8ce8abcf4304b4fbd Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Fri, 14 Jan 2011 13:49:36 +0100
Subject: [PATCH 5/9] Let PETScKrylovSolver report reason for divergence as string.
---
dolfin/la/PETScKrylovSolver.cpp | 5 +++--
1 files changed, 3 insertions(+), 2 deletions(-)
diff --git a/dolfin/la/PETScKrylovSolver.cpp b/dolfin/la/PETScKrylovSolver.cpp
index 373d396..ba3f3e6 100644
--- a/dolfin/la/PETScKrylovSolver.cpp
+++ b/dolfin/la/PETScKrylovSolver.cpp
@@ -201,11 +201,12 @@ dolfin::uint PETScKrylovSolver::solve(PETScVector& x, const PETScVector& b)
// Get solver residual norm
double rnorm = 0.0;
KSPGetResidualNorm(*_ksp, &rnorm);
+ const char *reason_str = KSPConvergedReasons[reason];
bool error_on_nonconvergence = parameters["error_on_nonconvergence"];
if (error_on_nonconvergence)
- error("PETSc Krylov solver did not converge in %i iterations (PETSc reason %i, norm %e).", num_iterations, reason, rnorm);
+ error("PETSc Krylov solver did not converge in %i iterations (PETSc reason %s, norm %e).", num_iterations, reason_str, rnorm);
else
- warning("Krylov solver did not converge in %i iterations (PETSc reason %i, norm %e).", num_iterations, reason, rnorm);
+ warning("Krylov solver did not converge in %i iterations (PETSc reason %s, norm %e).", num_iterations, reason_str, rnorm);
}
// Report results
--
1.7.1
From c3d524de2fc3cf4fa09a03db3d3c0791e3f2f284 Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Thu, 20 Jan 2011 10:36:36 +0100
Subject: [PATCH 6/9] Changes required to compile with petsc-dev.
---
cmake/modules/FindPETSc.cmake | 4 +++-
dolfin/la/PETScLUSolver.cpp | 12 ++++++++++++
dolfin/la/PETScMatrix.cpp | 18 ++++++++++++++++--
dolfin/la/PETScVector.cpp | 4 ++++
4 files changed, 35 insertions(+), 3 deletions(-)
diff --git a/cmake/modules/FindPETSc.cmake b/cmake/modules/FindPETSc.cmake
index 6875162..e226851 100644
--- a/cmake/modules/FindPETSc.cmake
+++ b/cmake/modules/FindPETSc.cmake
@@ -111,7 +111,9 @@ show :
endmacro()
# Call macro to get the PETSc variables
- petsc_get_variable(PETSC_INCLUDE PETSC_INCLUDE)
+ petsc_get_variable(PETSC_INCLUDE PETSC_INCLUDE) # 3.1
+ petsc_get_variable(PETSC_CC_INCLUDES PETSC_CC_INCLUDES) # dev
+ set(PETSC_INCLUDE ${PETSC_INCLUDE} ${PETSC_CC_INCLUDES})
petsc_get_variable(PETSC_LIB PETSC_LIB)
# Remove temporary Makefile
diff --git a/dolfin/la/PETScLUSolver.cpp b/dolfin/la/PETScLUSolver.cpp
index 1cd2afe..4264e4b 100644
--- a/dolfin/la/PETScLUSolver.cpp
+++ b/dolfin/la/PETScLUSolver.cpp
@@ -38,6 +38,18 @@ namespace dolfin
}
};
}
+
+// Forward compatibility with petsc-dev
+#if (PETSC_VERSION_RELEASE == 0)
+#define MAT_SOLVER_UMFPACK MATSOLVERUMFPACK
+#define MAT_SOLVER_MUMPS MATSOLVERMUMPS
+#define MAT_SOLVER_PASTIX MATSOLVERPASTIX
+#define MAT_SOLVER_PETSC MATSOLVERPETSC
+#define MAT_SOLVER_SPOOLES MATSOLVERSPOOLES
+#define MAT_SOLVER_SUPERLU_DIST MATSOLVERSUPERLU_DIST
+#define MAT_SOLVER_SUPERLU MATSOLVERSUPERLU
+#endif
+
//-----------------------------------------------------------------------------
// Available LU solver
const std::map<std::string, const MatSolverPackage> PETScLUSolver::lu_packages
diff --git a/dolfin/la/PETScMatrix.cpp b/dolfin/la/PETScMatrix.cpp
index f89ffba..3a219a0 100644
--- a/dolfin/la/PETScMatrix.cpp
+++ b/dolfin/la/PETScMatrix.cpp
@@ -288,10 +288,17 @@ void PETScMatrix::zero(uint m, const uint* rows)
assert(A);
IS is = 0;
+ PetscScalar null = 0.0;
+#if (PETSC_VERSION_RELEASE == 0)
+ ISCreateGeneral(PETSC_COMM_SELF, static_cast<int>(m),
+ reinterpret_cast<const int*>(rows),
+ PETSC_COPY_VALUES, &is);
+ MatZeroRowsIS(*A, is, null, NULL, NULL);
+#else
ISCreateGeneral(PETSC_COMM_SELF, static_cast<int>(m),
reinterpret_cast<const int*>(rows), &is);
- PetscScalar null = 0.0;
MatZeroRowsIS(*A, is, null);
+#endif
ISDestroy(is);
}
//-----------------------------------------------------------------------------
@@ -300,10 +307,17 @@ void PETScMatrix::ident(uint m, const uint* rows)
assert(A);
IS is = 0;
+ PetscScalar one = 1.0;
+#if (PETSC_VERSION_RELEASE == 0)
+ ISCreateGeneral(PETSC_COMM_SELF, static_cast<int>(m),
+ reinterpret_cast<const int*>(rows),
+ PETSC_COPY_VALUES, &is);
+ MatZeroRowsIS(*A, is, one, NULL, NULL);
+#else
ISCreateGeneral(PETSC_COMM_SELF, static_cast<int>(m),
reinterpret_cast<const int*>(rows), &is);
- PetscScalar one = 1.0;
MatZeroRowsIS(*A, is, one);
+#endif
ISDestroy(is);
}
//-----------------------------------------------------------------------------
diff --git a/dolfin/la/PETScVector.cpp b/dolfin/la/PETScVector.cpp
index 12a3df3..1b5de0d 100644
--- a/dolfin/la/PETScVector.cpp
+++ b/dolfin/la/PETScVector.cpp
@@ -612,7 +612,11 @@ void PETScVector::gather(GenericVector& y, const Array<uint>& indices) const
// Create local index sets
IS from, to;
+#if (PETSC_VERSION_RELEASE == 0)
+ ISCreateGeneral(PETSC_COMM_SELF, n, global_indices, PETSC_COPY_VALUES, &from);
+#else
ISCreateGeneral(PETSC_COMM_SELF, n, global_indices, &from);
+#endif
ISCreateStride(PETSC_COMM_SELF, n, 0 , 1, &to);
// Resize vector if required
--
1.7.1
From 9ac73fdb5c649e67b11952f78d2fa07d5fa0de2c Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Sun, 16 Jan 2011 14:13:07 +0100
Subject: [PATCH 7/9] NoDeleter doesn't require type information.
---
dolfin/common/Array.h | 4 ++--
dolfin/common/NoDeleter.h | 5 ++---
dolfin/la/PETScKrylovMatrix.cpp | 4 ++--
dolfin/la/PETScUserPreconditioner.cpp | 4 ++--
4 files changed, 8 insertions(+), 9 deletions(-)
diff --git a/dolfin/common/Array.h b/dolfin/common/Array.h
index af36738..cc72f5f 100644
--- a/dolfin/common/Array.h
+++ b/dolfin/common/Array.h
@@ -44,7 +44,7 @@ namespace dolfin
Array(uint N, boost::shared_array<T> x) : _size(N), x(x) {}
/// Construct array from a pointer. Array will not take ownership.
- Array(uint N, T* x) : _size(N), x(boost::shared_array<T>(x, NoDeleter<T>())) {}
+ Array(uint N, T* x) : _size(N), x(boost::shared_array<T>(x, NoDeleter())) {}
/// Destructor
~Array() {}
@@ -61,7 +61,7 @@ namespace dolfin
void update(uint N, T* _x)
{
_size = N;
- x.reset(_x, NoDeleter<T>());
+ x.reset(_x, NoDeleter());
}
/// Return informal string representation (pretty-print).
diff --git a/dolfin/common/NoDeleter.h b/dolfin/common/NoDeleter.h
index 098638d..ac520ae 100644
--- a/dolfin/common/NoDeleter.h
+++ b/dolfin/common/NoDeleter.h
@@ -16,11 +16,10 @@ namespace dolfin
/// NoDeleter is a customised deleter intended for use with smart pointers.
- template <class T>
class NoDeleter
{
public:
- void operator() (T* p) {}
+ void operator() (const void *p) {}
};
/// Helper function to construct shared pointer with NoDeleter with cleaner syntax
@@ -28,7 +27,7 @@ namespace dolfin
template<class T>
boost::shared_ptr<T> reference_to_no_delete_pointer(T& r)
{
- return boost::shared_ptr<T>(&r, NoDeleter<T>());
+ return boost::shared_ptr<T>(&r, NoDeleter());
}
}
diff --git a/dolfin/la/PETScKrylovMatrix.cpp b/dolfin/la/PETScKrylovMatrix.cpp
index 9c83349..c1d1dd9 100644
--- a/dolfin/la/PETScKrylovMatrix.cpp
+++ b/dolfin/la/PETScKrylovMatrix.cpp
@@ -28,8 +28,8 @@ namespace dolfin
int usermult(Mat A, Vec x, Vec y)
{
// Wrap x and y in a shared_ptr
- boost::shared_ptr<Vec> _x(&x, NoDeleter<Vec>());
- boost::shared_ptr<Vec> _y(&y, NoDeleter<Vec>());
+ boost::shared_ptr<Vec> _x(&x, NoDeleter());
+ boost::shared_ptr<Vec> _y(&y, NoDeleter());
void* ctx = 0;
MatShellGetContext(A, &ctx);
diff --git a/dolfin/la/PETScUserPreconditioner.cpp b/dolfin/la/PETScUserPreconditioner.cpp
index 82b1aea..450fb28 100644
--- a/dolfin/la/PETScUserPreconditioner.cpp
+++ b/dolfin/la/PETScUserPreconditioner.cpp
@@ -48,8 +48,8 @@ int PETScUserPreconditioner::PCApply(PC pc, Vec x, Vec y)
PETScUserPreconditioner* newpc = (PETScUserPreconditioner*)pc->data;
- boost::shared_ptr<Vec> _x(&x, NoDeleter<Vec>());
- boost::shared_ptr<Vec> _y(&y, NoDeleter<Vec>());
+ boost::shared_ptr<Vec> _x(&x, NoDeleter());
+ boost::shared_ptr<Vec> _y(&y, NoDeleter());
PETScVector dolfinx(_x), dolfiny(_y);
newpc->solve(dolfiny, dolfinx);
--
1.7.1
From 0d853a4f7b4aeb00dfcd4d93d89840511763ce0c Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Wed, 19 Jan 2011 22:03:25 +0100
Subject: [PATCH 8/9] Fix epetra crash with assembly of rectangular matrices.
---
dolfin/la/EpetraMatrix.cpp | 40 +++++++++++++++++++++++++++++++---------
dolfin/la/SparsityPattern.cpp | 2 +-
2 files changed, 32 insertions(+), 10 deletions(-)
diff --git a/dolfin/la/EpetraMatrix.cpp b/dolfin/la/EpetraMatrix.cpp
index fe34b55..1ccdeab 100644
--- a/dolfin/la/EpetraMatrix.cpp
+++ b/dolfin/la/EpetraMatrix.cpp
@@ -92,6 +92,14 @@ void EpetraMatrix::init(const GenericSparsityPattern& sparsity_pattern)
Epetra_MpiComm comm = f.get_mpi_comm();
Epetra_Map row_map(sparsity_pattern.size(0), num_local_rows, 0, comm);
+ // For rectangular matrices with more columns than rows, the columns which are
+ // larger than those in row_map are marked as nonlocal (and assembly fails).
+ // The domain_map fixes that problem, at least in the serial case.
+ // FIXME: Needs attention in the parallel case. Maybe range_map is also req'd.
+ const std::pair<uint,uint> colrange = sparsity_pattern.local_range(1);
+ const int num_local_cols = colrange.second - colrange.first;
+ Epetra_Map domain_map(sparsity_pattern.size(1), num_local_cols, 0, comm);
+
// Create Epetra_FECrsGraph
Epetra_FECrsGraph matrix_map(Copy, row_map, reinterpret_cast<int*>(&dnum_nonzeros[0]));
@@ -105,18 +113,32 @@ void EpetraMatrix::init(const GenericSparsityPattern& sparsity_pattern)
matrix_map.InsertGlobalIndices(global_row, row_set->size(), reinterpret_cast<int*>(&_nz_entries[0]));
}
- // Add off-diagonal block indices
- for (row_set = o_pattern.begin(); row_set != o_pattern.end(); ++row_set)
+ for (uint local_row = 0; local_row < d_pattern.size(); local_row++)
{
- const uint global_row = row_set - o_pattern.begin() + n0;
- const std::vector<dolfin::uint>& nz_entries = *row_set;
- std::vector<dolfin::uint>& _nz_entries = const_cast<std::vector<dolfin::uint>& >(nz_entries);
- matrix_map.InsertGlobalIndices(global_row, row_set->size(), reinterpret_cast<int*>(&_nz_entries[0]));
+ const uint global_row = local_row + n0;
+ std::vector<uint> &entries = const_cast<std::vector<uint>&>(d_pattern[local_row]);
+ matrix_map.InsertGlobalIndices(global_row, entries.size(), reinterpret_cast<int*>(&entries[0]));
}
- // Finalise map
- matrix_map.GlobalAssemble();
- matrix_map.OptimizeStorage();
+ // Add off-diagonal block indices (parallel only)
+ for (uint local_row = 0; local_row < o_pattern.size(); local_row++)
+ {
+ const uint global_row = local_row + n0;
+ std::vector<uint> &entries = const_cast<std::vector<uint>&>(o_pattern[local_row]);
+ matrix_map.InsertGlobalIndices(global_row, entries.size(), reinterpret_cast<int*>(&entries[0]));
+ }
+
+ try
+ {
+ // Finalise map. Here, row_map is standing in for RangeMap, which is
+ // probably ok but should be double-checked.
+ matrix_map.GlobalAssemble(domain_map, row_map);
+ matrix_map.OptimizeStorage();
+ }
+ catch (int err)
+ {
+ error("Epetra threw error %d in assemble", err);
+ }
// Create matrix
A.reset( new Epetra_FECrsMatrix(Copy, matrix_map) );
diff --git a/dolfin/la/SparsityPattern.cpp b/dolfin/la/SparsityPattern.cpp
index 0aee3ee..3af9e27 100644
--- a/dolfin/la/SparsityPattern.cpp
+++ b/dolfin/la/SparsityPattern.cpp
@@ -299,7 +299,7 @@ std::vector<std::vector<dolfin::uint> > SparsityPattern::diagonal_pattern(Type t
//-----------------------------------------------------------------------------
std::vector<std::vector<dolfin::uint> > SparsityPattern::off_diagonal_pattern(Type type) const
{
- std::vector<std::vector<uint> > v(diagonal.size());
+ std::vector<std::vector<uint> > v(off_diagonal.size());
for (uint i = 0; i < off_diagonal.size(); ++i)
v[i] = std::vector<uint>(off_diagonal[i].begin(), off_diagonal[i].end());
--
1.7.1
From fa98df7552078a59f5f6e2b52c5f928981c6c32b Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Thu, 20 Jan 2011 01:20:31 +0100
Subject: [PATCH 9/9] Implement DirichletBC::zero_columns for symmetric modification of off-diagonal block matrices.
---
dolfin/fem/DirichletBC.cpp | 54 ++++++++++++++++++++++++++++++++++++++++++++
dolfin/fem/DirichletBC.h | 4 +++
dolfin/la/EpetraMatrix.cpp | 10 +++++++-
3 files changed, 67 insertions(+), 1 deletions(-)
diff --git a/dolfin/fem/DirichletBC.cpp b/dolfin/fem/DirichletBC.cpp
index ad6dc16..75b7b17 100644
--- a/dolfin/fem/DirichletBC.cpp
+++ b/dolfin/fem/DirichletBC.cpp
@@ -209,6 +209,60 @@ void DirichletBC::zero(GenericMatrix& A) const
A.apply("insert");
}
//-----------------------------------------------------------------------------
+void DirichletBC::zero_columns(GenericMatrix& A, GenericVector& b) const
+{
+ Map bv_map;
+ get_boundary_values(bv_map, _method);
+
+ // Create lookup table of dofs
+ const uint nrows = A.size(0); // should be equal to b.size()
+ const uint ncols = A.size(1); // should be equal to max possible dof+1
+
+ std::vector<char> is_dof(ncols);
+ std::vector<double> dof_val(ncols);
+ for (Map::const_iterator bv = bv_map.begin(); bv != bv_map.end(); ++bv)
+ {
+ is_dof[bv->first] = 1;
+ dof_val[bv->first] = bv->second;
+ }
+
+ // Scan through all columns of all rows, setting to zero if is_dof[column]
+ // At the same time, we collect corrections to the RHS
+
+ std::vector<uint> cols;
+ std::vector<double> vals;
+ std::vector<double> b_vals;
+ std::vector<uint> b_rows;
+
+ for (uint row=0; row<nrows; row++)
+ {
+ A.getrow(row, cols, vals);
+ bool row_changed=false;
+ for (uint j=0; j<cols.size(); j++)
+ {
+ const uint col = cols[j];
+ if (is_dof[col] && vals[j] != 0.0)
+ {
+ if (!row_changed)
+ {
+ row_changed = true;
+ b_rows.push_back(row);
+ b_vals.push_back(0.0);
+ }
+ b_vals.back() -= dof_val[col]*vals[j];
+ vals[j] = 0.0;
+ }
+ }
+ if (row_changed) {
+ A.setrow(row, cols, vals);
+ A.apply("insert");
+ }
+ }
+
+ b.add(&b_vals.front(), b_rows.size(), &b_rows.front());
+ b.apply("add");
+}
+//-----------------------------------------------------------------------------
const std::vector<std::pair<dolfin::uint, dolfin::uint> >& DirichletBC::markers()
{
return facets;
diff --git a/dolfin/fem/DirichletBC.h b/dolfin/fem/DirichletBC.h
index 6cfcda8..cf15127 100644
--- a/dolfin/fem/DirichletBC.h
+++ b/dolfin/fem/DirichletBC.h
@@ -169,6 +169,10 @@ namespace dolfin
/// non-diagonal matrices in a block matrix.
void zero(GenericMatrix& A) const;
+ /// Make columns associated with boundary conditions zero, and
+ /// update the RHS to reflect the changes. Useful for non-diagonals.
+ void zero_columns(GenericMatrix& A, GenericVector& b) const;
+
/// Return boundary markers (facets stored as pairs of cells and local
/// facet numbers)
const std::vector<std::pair<uint, uint> >& markers();
diff --git a/dolfin/la/EpetraMatrix.cpp b/dolfin/la/EpetraMatrix.cpp
index 1ccdeab..dfc2f56 100644
--- a/dolfin/la/EpetraMatrix.cpp
+++ b/dolfin/la/EpetraMatrix.cpp
@@ -491,7 +491,15 @@ void EpetraMatrix::getrow(uint row, std::vector<uint>& columns,
void EpetraMatrix::setrow(uint row, const std::vector<uint>& columns,
const std::vector<double>& values)
{
- dolfin_not_implemented();
+ static bool print_msg_once=true;
+ if (print_msg_once)
+ {
+ info("EpetraMatrix::setrow is implemented inefficiently");
+ print_msg_once = false;
+ }
+
+ for (uint i=0; i<columns.size(); i++)
+ set(&values[i], 1, &row, 1, &columns[i]);
}
//-----------------------------------------------------------------------------
LinearAlgebraFactory& EpetraMatrix::factory() const
--
1.7.1
Follow ups
