initial import

pyfefp/Brinkman.py

+import numpy as np
+from . import quadrature
+
+_u_ = 0
+_v_ = 1
+_w_ = 2
+
+def derivate(f, x, nd = 1, eps = 1e-8) :
+    f0 = f(x)
+    xd = x.reshape(x.shape[:x.ndim-nd]+(-1,))
+    d = np.ndarray((f0.shape +(xd.shape[-1],)))
+    for i in range(xd.shape[-1]) :
+        xd[..., i] += eps
+        d[..., i] = (f(xd.reshape(x.shape)) - f0) / eps
+        xd[..., i] -= eps
+    return d.reshape(f0.shape + x.shape[x.ndim - nd:])
+
+
+class Brinkman :
+
+    def dot(self, a, b) :
+        return a[..., 0] * b[..., 0] + a[..., 1] * b[..., 1] + (0 if self.D == 2 else a[..., 2] * b[..., 2])
+
+    def psiTerm0(self, sol, solgrad):
+        _p_ = self.D
+        _U_ = range(self.D)
+        val = np.zeros_like(sol)
+        #self.g[None, None, :] * self.ca * self.rhof
+        val[..., _U_] =   - (sol[...,_U_] - self.sol0qp[...,_U_]) * self.rhof/self.dt\
+            - solgrad[..., _p_, :] \
+            + self.mu * (self.dot(solgrad[...,_U_,:], self.dca)-sol[...,_U_]/self.ca**2*self.dot(self.dca, self.dca))
+        val[..., _p_] = - solgrad[..., _u_, 0] - solgrad[..., _v_, 1] - (0 if self.D == 2 else solgrad[..., _w_, 2])
+        val[..., _p_] += - self.dvs[..., 0, 0] - self.dvs[..., 1, 1] - (0 if self.D == 2 else self.dvs[..., 2, 2])
+        #val[..., _v_] += self.xqp[...,0] * (1-self.xqp[...,0]**2/4)
+        return val
+
+    def psiTerm1(self, sol, solgrad):
+        _p_ = self.D
+        _U_ = range(self.D)
+
+        val = np.ndarray(sol.shape + (self.D,))
+        val[..., _U_, :] = -solgrad[..., _U_, :] * self.mu
+        val[..., _p_, :] =  0#self.vs# + sol[..., _U_]
+        val[..., _p_, :] = -1e-8 * solgrad[..., _p_, :]
+        if self.ns :
+            val[..., _u_, :] += - sol[..., [_u_]] * sol[..., _U_] / self.ca * self.rhof
+            val[..., _v_, :] += - sol[..., [_v_]] * sol[..., _U_] / self.ca * self.rhof
+            if self.D == 3:
+                val[..., _w_, :] += - sol[..., [_w_]] * sol[..., _U_] / self.ca * self.rhof
+        return val
+
+    # forces by the particle ON the fluid
+    def particleForces(self, sol, solgrad, fonpart = None) :
+        _p_ = self.D
+        _U_ = range(self.D)
+        val = np.zeros_like(sol)
+        if self.D == 2 :
+            d = 2 * (self.pf.volume[0, 0] / np.pi)**0.5
+        else :
+            d = 2 * (3./(4 * np.pi) * self.pf.volume[0, 0])**(1./3)
+        U = self.vs - sol[..., _U_]/self.ca
+        normU = np.sqrt(self.dot(U, U))[...,np.newaxis]
+        Re_O_u = d * self.ca/self.mu * self.rhof
+        Cd_u = (0.63 * normU + 4.8/(Re_O_u**0.5))**2
+        f = self.ca**(-1.8)
+        F = solgrad[..., _p_, :]
+        GoU = f * Cd_u * self.rhof / 2 * d
+        #mass = self.pf.volume * 3/2 * self.rhop
+        aext = self.g * (self.rhop - self.rhof) /self.rhop
+        UImp = U + aext * self.stab
+        if fonpart is not None:
+            fonpart[...] = (-F*self.pf.volume - GoU*UImp) / (1 + self.stab/self.pf.mass * GoU) + aext * self.pf.mass
+        val[..., _U_] =  (F*self.pf.volume + GoU* UImp) / (1 + self.stab/self.pf.mass * GoU)
+        val[..., _p_] = 0
+
+        return val
+
+    def __init__(self, mu, pf, g, rho, rhop, ns = False, stab = 0, D = 2) :
+        self.D = D
+        self.ns = ns
+        self.rhof = rho
+        self.rhop = rhop
+        self.stab = stab
+        self.mu = mu
+        self.pf = pf
+        self.g = np.array(g)
+
+    def particleInteractionTerm(self, meshJ, solution, res, jac) :
+        _p_ = self.D
+        _U_ = range(self.D)
+        dofm = solution.dof_manager()
+        self.ca = self.pf.porosity.at_points(self.pf.eid, self.pf.uvw)
+        self.vs = self.pf.velocity
+        sol = solution.at_points(self.pf.eid, self.pf.uvw)
+        gradsol = solution.gradient_at_points(self.pf.eid, self.pf.uvw, meshJ)
+        pTerm = self.particleForces(sol, gradsol, self.pf.forces)
+        for i in range(dofm.n_fields()) :
+            np.add.at(res, dofm._fields[i][0][:, self.pf.eid], (dofm.get_field_basis(i).f(self.pf.uvw)*pTerm[:, [i]]).T)
+        d00 = derivate(lambda s : self.particleForces(s, gradsol), sol, 1)
+        d01 = derivate(lambda s : self.particleForces(sol, s), gradsol, 2)
+        for f in range(dofm.n_fields()) :
+            ff = dofm.get_field_basis(f).f(self.pf.uvw)
+            for g in range(dofm.n_fields()) :
+                pp = d00[:, f, g, None] * dofm.get_field_basis(g).f(self.pf.uvw)
+                d01xi = meshJ.multDXiDXLeft(d01[:, f, g, :], 1, self.pf.eid)
+                pp += (d01xi[:, None, :] * dofm.get_field_basis(g).df(self.pf.uvw)).sum(2)
+                jac.add_to_field(f, g, pp[:, None, :] * ff[:, :, None], self.pf.eid)
+
+    def fluidTerm(self, meshJ, solution, res, jac) :
+        dofm = solution.dof_manager()
+        qp, qw = quadrature.get_triangle(3) if self.D == 2 else quadrature.get_tetrahedron(3)
+        self.vs = self.pf.momentum.at_qp(qp)
+        self.dvs = self.pf.momentum.gradient_at_qp(qp, meshJ)
+        self.ca = self.pf.porosity.at_qp(qp)
+        self.dca = self.pf.porosity.gradient_at_qp(qp, meshJ)
+        sol = solution.at_qp(qp)
+        gradsol = solution.gradient_at_qp(qp, meshJ)
+        self.sol0qp = self.sol0.at_qp(qp)
+        #f = meshJ.dofManager.get_field_basis(0).f(qp)
+        #self.xqp = np.tensordot(self.x,f,[1,1]).swapaxes(1,2)
+
+        x0 = meshJ.multJ(self.psiTerm0(sol, gradsol))
+        x1 = meshJ.multJ(meshJ.multDXiDXLeft(self.psiTerm1(sol, gradsol), 3))
+        for i in range(dofm.n_fields()) :
+            basis = dofm.get_field_basis(i)
+            psiW = basis.f(qp) * qw[:, None]
+            r0 = np.tensordot(psiW, x0[:, :, i], [0, 1])
+            dPsiW = basis.df(qp) * qw[:, None, None]
+            r1 = np.tensordot(dPsiW, x1[:, :, i, :], [[0, 2], [1, 2]])
+            np.add.at(res, dofm._fields[i][0], r0 + r1)
+
+        d00 = derivate(lambda s : self.psiTerm0(s, gradsol), sol, 1)
+        d01 = derivate(lambda s : self.psiTerm0(sol, s), gradsol, 2)
+        d10 = derivate(lambda s : self.psiTerm1(s, gradsol) , sol, 1)
+        d11 = derivate(lambda s : self.psiTerm1(sol, s), gradsol, 2)
+        for f in range(dofm.n_fields()):
+            basisf = dofm.get_field_basis(f)
+            for g in range(dofm.n_fields()):
+                basisg = dofm.get_field_basis(g)
+                if (np.max(np.abs(d00[...,f,g])) > 1e-8) :
+                    psiPsiW = np.einsum("qn, qm, q-> qnm", basisf.f(qp), basisg.f(qp), qw)
+                    jac.add_to_field(f, g, np.tensordot(meshJ.multJ(d00[:,:,f,g]), psiPsiW, 1))
+                if (np.max(np.abs(d01[:,:,f,g,:])) > 1e-8) :
+                    dPsiPsiW = np.einsum("qmd, qn, q -> qdnm", basisg.df(qp), basisf.f(qp), qw)
+                    jac.add_to_field(f, g, np.tensordot(meshJ.multJ(meshJ.multDXiDXLeft(d01[:,:,f,g,:], 2)), dPsiPsiW, 2))
+                if (np.max(np.abs(d10[:,:,f,:,g])) > 1e-8) :
+                    psiDPsiW = np.einsum("qmd, qn, q -> qdmn", basisf.df(qp), basisg.f(qp), qw)
+                    jac.add_to_field(f, g, np.tensordot(meshJ.multJ(meshJ.multDXiDXLeft(d10[:,:,f,:,g], 2)), psiDPsiW, 2))
+                if (np.max(np.abs(d11[:,:,f,:,g,:])) > 1e-8) :
+                    dPsiDPsiW = np.einsum("qmx, qny, q -> qxymn", basisf.df(qp), basisg.df(qp), qw)
+                    jac.add_to_field(f, g, np.tensordot(meshJ.multDXiDXLeft(meshJ.multDXiDXLeft(meshJ.multJ(d11[:,:,f,:,g,:]), 3), 2), dPsiDPsiW, 3))
+
+    def massTerm(self, alpha, solution, res, jac) :
+        pass
+
+    def volumeTerm(self, meshJ, solution) :
+        dofm = solution.dof_manager()
+        jac = dofm.new_matrix()
+        res = dofm.new_vector()
+        self.particleInteractionTerm(meshJ, solution, res, jac)
+        self.fluidTerm(meshJ, solution, res, jac)
+        #self.massTerm(meshJ, alphamass, solution, res, jac)
+        return res, jac.to_csr_matrix()

pyfefp/Steady.py

+import numpy as np
+import time
+import scipy.sparse
+import scipy.sparse.linalg
+
+class Steady :
+    def __init__(self, law, dof, jac, ilu=True) :
+        self.law = law
+        self.meshJacobian = jac
+        self.dof = dof
+        self.preconditioner = None
+        self.tspsolve = 0
+        self.tvt = 0
+        self.tas = 0
+        self.iterative = True
+        self.ilu = ilu
+        self.x = self.dof.new_vector()
+
+    def solve(self, rtol = 1e-8, atol = 1e-10) :
+        x = self.x
+        self.dof.update_boundary_conditions(self.meshJacobian)
+        x[self.dof._fixedidx] = self.dof._fixedval
+        firstR = None
+        newtonIter = 0
+        ttot = time.clock()
+        while True :
+            tic = time.clock()
+            res, jac = self.law.volumeTerm(self.meshJacobian, x)
+            self.tvt += time.clock() - tic
+
+            tic = time.clock()
+            for row in self.dof._fixedidx :
+                a = slice(jac.indptr[row], jac.indptr[row +1])
+                jac.data[a] = np.where(jac.indices[a] == row, 1., 0.)
+            jac.eliminate_zeros()
+            self.tas += time.clock() - tic
+            res[self.dof._fixedidx] = 0
+            resn = np.linalg.norm(res)
+            if firstR :
+                print("newton (%i) % 6.2e %.1f" % (newtonIter, resn, -np.log10(resn/firstR)))
+                if (resn/firstR < rtol or resn < atol) :
+                    break
+            else :
+                firstR = resn
+                print("newton (%i) % 6.2e" % (newtonIter, resn))
+            newtonIter += 1
+            if newtonIter > 10:
+                break
+            tic = time.clock()
+            matrix = jac.tocsc()
+            if not self.preconditioner :
+                P = scipy.sparse.linalg.splu(matrix)
+                self.preconditioner = scipy.sparse.linalg.LinearOperator(matrix.shape, lambda x: P.solve(x))
+            if self.iterative :
+                info = 1
+                if self.preconditioner :
+                    #I think there is a bug in scipy, without the callback argument, maxiter as no effect
+                    dx, info = scipy.sparse.linalg.gmres(matrix, res, tol=1e-10, M=self.preconditioner, maxiter=50, callback=lambda x:None)
+                #if info != 0 and self.ilu:
+                #    print("*** PRECOND ***")
+                #    P = scipy.sparse.linalg.spilu(matrix)
+                #    self.preconditioner = scipy.sparse.linalg.LinearOperator(matrix.shape, lambda x: P.solve(x))
+                #    dx,info  = scipy.sparse.linalg.gmres(matrix, res, tol=1e-10, M=self.preconditioner, maxiter=50, callback=lambda x:None)
+                if info != 0 :
+                    print("*** PRECOND FULL LU***")
+                    P = scipy.sparse.linalg.splu(matrix)
+                    self.preconditioner = scipy.sparse.linalg.LinearOperator(matrix.shape, lambda x: P.solve(x))
+                    dx,info  = scipy.sparse.linalg.gmres(matrix, res, tol=1e-10, M=self.preconditioner, maxiter=15, callback=lambda x:None)
+            else :
+                dx = scipy.sparse.linalg.spsolve(matrix, res)
+            x -= dx
+            self.tspsolve += time.clock() - tic
+
+        print("spsolve: %.3g  term: %.3g  assemble: %.3g" % (self.tspsolve, self.tvt, self.tas))
+        return x

pyfefp/__init__.py

+from .Brinkman import Brinkman
+from .Steady import Steady
+from .dofManager import DofManager
+from .meshJacobian import MeshJacobian
+from .legendre import *
+from .particleFluid import ParticleFluid
+from .fluidSolver import FluidSolver
+from .mesh import Mesh

pyfefp/dofManager.py

+from . import mesh
+import numpy as np
+import scipy.sparse
+
+class DofManager :
+
+    def get_n(self) :
+        return self._ndof
+
+    def get_field_basis(self, i) :
+        return self._fields[i][2]
+
+    def getField(self, i, v) :
+        return v[self._fields[i][0]]
+
+    def getElements(self, iField) :
+        return self._fields[iField][1]
+
+    def addToField(self, i, vf, v) :
+        np.add.at(v, self._fields[i][0], vf)
+
+    def setField(self, i, vf, v) :
+        v[self._fields[i][0]] =  vf
+
+    def n_fields(self) :
+        return len(self._fields)
+
+    def add_field(self, basis) :
+        field = []
+        elist = []
+        ifield = len(self._fields)
+        vdof = self._dof_by_dimension[0]
+        nV = basis.n(0)
+        for v in self._mesh.vertices :
+            vdof[(ifield, v[3])] = vdof.get((ifield, v[3]), ()) + tuple(range(self._ndof, self._ndof + nV))
+            self._ndof += nV
+        for en in self._mesh.entities :
+            n = basis.n(en.dimension)
+            dof = self._dof_by_dimension[en.dimension]
+            if en.dimension != basis.dimension() :
+                continue
+            for e in en.elements :
+                dof[(ifield, e.tag)] = dof.get((ifield, e.tag), ()) + tuple(range(self._ndof, self._ndof + n))
+                self._ndof += n
+                elist.append(e)
+                cdof = ()
+                for v in e.vertices :
+                    cdof += vdof.get((ifield, v[3]), ())
+                cdof += dof[(ifield, e.tag)]
+                field.append(cdof)
+        self._fields.append((np.transpose(np.array(field, np.int)), elist, basis))
+
+    def add_boundary_condition(self, tag, field, val = None, cb = None) :
+        if cb is not None :
+            val = 0.
+        vdof = self._dof_by_dimension[0]
+        fixed_start = len(self._fixedidx)
+        fixedvid = []
+        for en in self._mesh.entities :
+            for ph in en.physicals :
+                phname = self._mesh.getPhysicalName(en.dimension, ph)
+                if phname != tag :
+                    continue
+                for e in en.elements :
+                    for x, y, z, vid in e.vertices :
+                        for dof in vdof.get((field, vid), ()) :
+                            if not dof in self._fixed :
+                                self._fixed[dof] = len(self._fixedval)
+                                self._fixedval.append(val)
+                                self._fixedidx.append(dof)
+                                fixedvid.append(vid)
+        if cb is not None :
+            self._boundary_cb.append((cb, slice(fixed_start, len(self._fixedidx)), fixedvid))
+
+
+    class matrix(np.ndarray) :
+
+        def __new__(cls, csrrows, slicedidx) :
+            obj = np.zeros(csrrows.shape).view(cls)
+            obj._csrrows = csrrows
+            obj._slicedidx = slicedidx
+            return obj
+
+        def to_csr_matrix(self) :
+            return scipy.sparse.csr_matrix((self.flat, (self._csrrows.flat, self._csrrows.swapaxes(1, 2).flat)))
+
+        def add_to_field(self, f, g, v, els = None) :
+            if els is None:
+                self[:, self._slicedidx[f]:self._slicedidx[f+1], self._slicedidx[g]:self._slicedidx[g+1]] += v
+            else :
+                np.add.at(self[:, self._slicedidx[f]:self._slicedidx[f+1], self._slicedidx[g]:self._slicedidx[g+1]], els, v)
+
+
+    class vector(np.ndarray) :
+
+        def __new__(cls, dofmanager) :
+            obj = np.zeros(dofmanager.get_n()).view(cls)
+            obj._dofmanager = dofmanager
+            return obj
+
+        def dof_manager(self) :
+            return self._dofmanager
+
+        def at_qp(self, q) :
+            dofmanager = self._dofmanager
+            NE = len(dofmanager._fields[0][1])
+            val = np.ndarray((NE, q.shape[0], dofmanager.n_fields()))
+            for i in range(dofmanager.n_fields()) :
+                basis = dofmanager.get_field_basis(i)
+                val[..., i] = np.tensordot(dofmanager.getField(i, self), basis.f(q), [0, 1])
+            return val
+
+        def gradient_at_qp(self, q, j) :
+            dofmanager = self._dofmanager
+            NE = len(dofmanager._fields[0][1])
+            val = np.ndarray((NE, q.shape[0], dofmanager.n_fields(), j.dim))
+            for i in range(dofmanager.n_fields()) :
+                basis = dofmanager.get_field_basis(i)
+                valxi = np.tensordot(dofmanager.getField(i, self), basis.df(q) , [0, 1])
+                val[..., i, :] = j.multDXiDXRight(valxi, 2)
+            return val
+
+        def at_points(self, els, xi) :
+            dofmanager = self._dofmanager
+            val = np.ndarray((len(els), dofmanager.n_fields()))
+            for i in range(dofmanager.n_fields()) :
+                f = dofmanager.getField(i, self)[:, els]
+                basis = dofmanager.get_field_basis(i)
+                val[:, i] =  (f * basis.f(xi).T).sum(0)
+            return val
+
+        def gradient_at_points(self, els, xi, j) :
+            dofmanager = self._dofmanager
+            val = np.ndarray((len(els), dofmanager.n_fields(), j.dim))
+            for i in range(dofmanager.n_fields()) :
+                f = dofmanager.getField(i, self)[:, els]
+                basis = dofmanager.get_field_basis(i)
+                dxi = (basis.df(xi) * f.T[:, :, None]).sum(1)
+                val[:, i, :] = j.multDXiDXRight(dxi, 1, els)
+            return val
+
+
+    def new_vector(self) :
+        return self.vector(self)
+
+    def new_matrix(self) :
+        return self.matrix(self._jaccsrrows, self._slicedidx)
+
+    def update_boundary_conditions(self, jac) :
+        x = [jac.dofManager.new_vector() for i in range(jac.dim)]
+        for i in range(jac.dim) :
+            jac.dofManager.setField(0, jac.x[:, :, i].T, x[i])
+        for cb, idx, vid in self._boundary_cb :
+            vidx = [jac.dofManager._dof_by_dimension[0][(0, i)][0] for i in vid]
+            self._fixedval[idx] = cb(np.hstack((x[i][vidx, None] for i in range(jac.dim))))
+
+    def complete(self) :
+        self._fixedidx = np.array(self._fixedidx)
+        self._fixedval = np.array(self._fixedval)
+        self._slicedidx = [0] * (self.n_fields() + 1)
+        n_tot = sum((self.get_field_basis(i).n() for i in range(self.n_fields())))
+        self._jaccsrrows = np.ndarray((len(self.getElements(0)), n_tot, n_tot))
+        for f in range(self.n_fields()):
+            self._slicedidx[f + 1] = self._slicedidx[f] + self.get_field_basis(f).n()
+            self._jaccsrrows[:, self._slicedidx[f]:self._slicedidx[f + 1], :] = self._fields[f][0].T[..., None]
+
+    def __init__(self, mesh) :
+        self._dof_by_dimension = [{},{},{},{}]
+        self._mesh = mesh
+        self._ndof = 0
+        self._fields = []
+        self._fixedval = []
+        self._fixedidx = []
+        self._fixed = {}
+        self._boundary_cb = []

pyfefp/fluidSolver.py

+from .meshJacobian import MeshJacobian
+from .particleFluid import ParticleFluid
+from .Brinkman import Brinkman
+from .dofManager import DofManager
+from .Steady import Steady
+from .mesh import Mesh
+from .legendre import *
+import time
+
+class FluidSolver :
+    def __init__(self, mesh, dt, rhop, rho, mu, g, imex=True) :
+        if isinstance(mesh, str):
+            mesh = Mesh(mesh)
+        self._mesh = mesh
+        self._jac = MeshJacobian(mesh)
+        self._pf = ParticleFluid(self._jac)
+        self._imex = imex
+        self._dt = dt
+        self._stabFactor = 1
+        law = Brinkman(
+            mu = mu,
+            pf = self._pf,
+            ns = False,
+            g = np.array(g),
+            rho = rho,
+            rhop = rhop,
+            D = self._jac.dim,
+            stab = dt)
+        d = DofManager(mesh)
+        if self._jac.dim == 2 :
+            d.add_field(TriangleP1)
+            d.add_field(TriangleP1)
+            d.add_field(TriangleP1)
+        elif self._jac.dim == 3:
+            d.add_field(TetrahedronP1Mini)
+            d.add_field(TetrahedronP1Mini)
+            d.add_field(TetrahedronP1Mini)
+            d.add_field(TetrahedronP1)
+        else :
+            raise ValueError("invalid mesh dimension : %i" % self._jac.dim)
+        self._dt = dt
+        self._law = law
+        law.dt = dt
+        law.x = self._jac.x
+        self._dof = d
+        self._sol = None
+
+    def add_boundary_condition(self, tag, field, *a, **b) :
+        self._dof.add_boundary_condition(tag, field, *a, **b)
+
+    def complete(self) :
+        self._dof.complete()
+        self._sol = self._dof.new_vector()
+        self._law.sol0 = self._dof.new_vector()
+        self._sol[:] = 0
+        self._law.sol0[:] = 0
+        self._solver = Steady(self._law, self._dof, self._jac)
+        self._solverEx = Steady(self._law, self._dof, self._jac)
+
+    def solve(self) :
+        tic = time.clock()
+        if self._imex :
+            v0 = self._pf.velocity.copy()
+            self._law.stab = self._dt *0.5 * self._stabFactor
+            self._sol[:] = self._solver.solve()
+            self._pf.velocity += self._pf.forces * self._dt*0.5 / self._pf.mass
+            self._law.stab = 0
+            self._sol[:] = self._solverEx.solve()
+            self._pf.velocity[:] = v0
+        else :
+            self._law.stab = self._dt * self._stabFactor
+            self._sol[:] = self._solver.solve()
+        self._law.sol0[:] = self._sol
+
+        print("cpu : ", time.clock() - tic)
+        return self.forces
+
+    def set_stab_factor(self, s) :
+        self._stabFactor = s
+
+    def write_position(self, odir, oiter, time, filetype="msh"):
+        if filetype == "msh" :
+            basename = "%s/f_%05i" % (odir, oiter)
+            if self._jac.dim == 2 :
+                X = np.zeros(((self._jac.dim + 1) *3, self._jac.x.shape[0]))
+                X[[0,1,3,4,6,7], :] = self._jac.x.reshape([-1, 6]).T
+            else :
+                X = self._jac.x.reshape(-1, 12).T
+            self._jac.writeScalar(basename + "-xmesh.msh", "x", oiter, time, X)
+        else :
+            print("unknown export format : \"%s\"." % filetype)
+
+    def write_solution(self, odir, oiter, time, filetype="msh") :
+        dim = self._jac.dim
+        if filetype == "msh" :
+            basename = "%s/f_%05i" % (odir, oiter)
+            self._jac.writeScalar(basename + "-porosity.msh", "porosity", oiter, time, self._jac.dofManager.getField(0, self._pf.porosity))
+            self._jac.writeScalar(basename + "-pressure.msh", "pressure", oiter, time, self._dof.getField(dim, self._sol))
+            X = np.zeros(((dim+1)*3, self._jac.x.shape[0]))
+            for d in range (dim) :
+                X[d::3, :] = self._dof.getField(d, self._sol)[:-1,:]
+            self._jac.writeScalar(basename + "-velocity.msh", "velocity", oiter, time, X)
+        elif filetype == "vtk" :
+            output = open("%s/f_%05i.vtu" %(odir, oiter), "w")
+            vertices = np.array(self._mesh.vertices)[:,:3]
+            elements = np.array([[v[3] for v in e.vertices] for e in self._dof.getElements(0)], np.int32) -1
+            v = np.zeros((vertices.shape[0], 3))
+            for d in range (dim) :
+                v[elements, d] = self._jac.x[:,:,d]
+            output.write("<VTKFile type=\"UnstructuredGrid\" format=\"0.1\">\n")
+            output.write("<UnstructuredGrid>\n")
+            output.write("<Piece NumberOfPoints=\"%i\" NumberOfCells=\"%i\">\n" % (vertices.shape[0], elements.shape[0]))
+            output.write("<Points>\n")
+            output.write("<DataArray NumberOfComponents=\"3\" type=\"Float64\" format=\"ascii\">\n")
+            np.savetxt(output, v)
+            output.write("</DataArray>\n")
+            output.write("</Points>\n")
+            output.write("<Cells>\n")
+            output.write("<DataArray  type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n")
+            np.savetxt(output, elements, fmt="%i")
+            output.write("</DataArray>\n")
+            output.write("<DataArray  type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n")
+            np.savetxt(output, (np.array(range(elements.shape[0]), np.int32)+1)*(dim+1), fmt="%i", newline=" ")
+            output.write("</DataArray>\n")
+            output.write("<DataArray  type=\"Int32\" Name=\"types\" format=\"ascii\">\n")
+            t = np.ndarray((elements.shape[0]), np.int32)
+            t[:] = 5 if dim == 2 else 10
+            np.savetxt(output, t, fmt="%i", newline=" ")
+            output.write("</DataArray>\n")
+            output.write("</Cells>\n")
+            output.write("<PointData Scalars=\"Pressure Porosity\" Vectors=\"Velocity\">\n")
+            output.write("<DataArray Name=\"Porosity\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">")
+            np.savetxt(output, self._pf.porosity)
+            output.write("</DataArray>\n")
+            output.write("<DataArray Name=\"Pressure\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">")
+            p = np.ndarray((vertices.shape[0]))
+            p[np.transpose(elements)] = self._dof.getField(dim, self._sol)
+            np.savetxt(output, p)