Commit e01591ec authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

3d + update all testcases

parent 2044932b
import numpy as np
from . import quadrature
D = 2
_u_ = 0
_v_ = 1
_w_ = 2
_p_ = D
_U_ = range(D)
def derivate(f, x, nd = 1, eps = 1e-8) :
f0 = f(x)
......@@ -18,79 +15,95 @@ def derivate(f, x, nd = 1, eps = 1e-8) :
xd[..., i] -= eps
return d.reshape(f0.shape + x.shape[x.ndim - nd:])
def dot(a, b) :
return a[..., 0] * b[..., 0] + a[..., 1] * b[..., 1] + (0 if D == 2 else a[..., 2] * b[..., 2])
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)
val[..., _U_] = - solgrad[..., _p_, :] * self.ca + self.mu * (\
+ dot(solgrad[..., _U_, :], self.dca) \
- sol[..., _U_] / self.ca **2 * dot(self.dca, self.dca))
val[..., [_v_]] += -self.ca
val[..., _p_] = - solgrad[..., _u_, 0] - solgrad[..., _v_, 1] - (0 if D == 2 else solgrad[..., _w_, 2])
val[..., _U_] = \
self.g[None, None, :] * self.ca * self.rhof\
- 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])
return val
def psiTerm1(self, sol, solgrad):
val = np.ndarray(sol.shape + (D,))
_p_ = self.D
_U_ = range(self.D)
val = np.ndarray(sol.shape + (self.D,))
val[..., _U_, :] = -solgrad[..., _U_, :] * self.mu
val[..., _p_, :] = self.vs
val[..., _p_, :] = 0#self.vs# + sol[..., _U_]
if self.stab != 0 :
pforces = solgrad[..., _p_, :]#self.particleForces(sol, solgrad)[..., _U_]
val[..., _p_, :] -= pforces * self.stab * (1-self.ca)
if self.ns :
val[..., _u_, :] += - sol[..., [_u_]] * sol[..., _U_] / self.ca
val[..., _v_, :] += - sol[..., [_v_]] * sol[..., _U_] / self.ca
if D == 3:
val[..., _w_, :] += - sol[..., [_w_]] * sol[..., _U_] / self.ca
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
def particleTerm(self, sol, solgrad) :
# forces by the particle ON the fluid
def particleForces(self, sol, solgrad) :
_p_ = self.D
_U_ = range(self.D)
val = np.zeros_like(sol)
if D == 2 :
d = 2 * (self.pf.volume[0, 0] / np.pi)**0.5
if self.D == 2 :
d = 2 * (self.pf.volume[0, 0] * 3./2 / 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(dot(U, U))[...,np.newaxis]
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_, :]
G = f * Cd_u * self.rhof * U / 2
val[..., _U_] = G
val[..., _U_] = F + G
val[..., _p_] = 0
return val
def __init__(self, mu, pf, ns = False, rhof = 1) :
def __init__(self, mu, pf, g, ns = False, rho = 1, stab = 0, D = 2) :
self.D = D
self.ns = ns
self.rhof = 1
self.rhof = rho
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.particleTerm(sol, gradsol) * self.pf.volume
self.pf.forces[:, :] = -pTerm[:, _U_] - gradsol[..., _p_, :] * self.pf.volume - self.pf.volume * np.array([0, 1])
pTerm = self.particleForces(sol, gradsol) * self.pf.volume
self.pf.forces[:, :] = -(pTerm[:, _U_])
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.particleTerm(s, gradsol) * self.pf.volume, sol, 1)
d01 = derivate(lambda s : self.particleTerm(sol, s) * self.pf.volume, gradsol, 2)
d00 = derivate(lambda s : self.particleForces(s, gradsol), sol, 1) * self.pf.volume[..., None]
d01 = derivate(lambda s : self.particleForces(sol, s), gradsol, 2) * self.pf.volume[..., None, None]
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)
if (np.max(np.abs(d01[:,f,g,:])) > 1e-8) :
d01xi = meshJ.multDXiDXLeft(d01[:, f, g, :], 1, self.pf.eid)
pp += (d01xi[:, None, :] * dofm.get_field_basis(g).df(self.pf.uvw)).sum(2)
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)
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)
......@@ -129,10 +142,14 @@ class Brinkman :
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()
......@@ -4,7 +4,7 @@ import scipy.sparse
import scipy.sparse.linalg
class Steady :
def __init__(self, law, dof, jac) :
def __init__(self, law, dof, jac, ilu=True) :
self.law = law
self.meshJacobian = jac
self.dof = dof
......@@ -12,9 +12,13 @@ class Steady :
self.tspsolve = 0
self.tvt = 0
self.tas = 0
self.iterative = True
self.ilu = ilu
self.x = self.dof.new_vector()
def solve(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
......@@ -34,7 +38,7 @@ class Steady :
resn = np.linalg.norm(res)
if firstR :
print("newton (%i) % 6.2e %.1f" % (newtonIter, resn, -np.log10(resn/firstR)))
if (resn/firstR < 1e-6) :
if (resn/firstR < rtol or resn < atol) :
break
else :
firstR = resn
......@@ -45,15 +49,25 @@ class Steady :
tic = time.clock()
matrix = jac.tocsc()
if not self.preconditioner :
P = scipy.sparse.linalg.spilu(matrix)
P = scipy.sparse.linalg.splu(matrix)
self.preconditioner = scipy.sparse.linalg.LinearOperator(matrix.shape, lambda x: P.solve(x))
#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=15, callback=lambda x:None)
if info != 0 :
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)
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
......
......@@ -4,12 +4,9 @@ import scipy.sparse
class DofManager :
def getN(self) :
def get_n(self) :
return self._ndof
def getNFree(self) :
return self._ndof - len(self._fixedidx)
def get_field_basis(self, i) :
return self._fields[i][2]
......@@ -53,25 +50,12 @@ class DofManager :
field.append(cdof)
self._fields.append((np.transpose(np.array(field, np.int)), elist, basis))
def set_boundary_condition(self, tag, field, val) :
vdof = self._dof_by_dimension[0]
edof = self._dof_by_dimension[self._fields[field][2].dimension()]
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), ()) :
self._fixedval[self._fixed[dof]] = val
for dof in edof.get((field, e.tag), ()) :
self._fixedval[self._fixed[dof]] = val
def add_boundary_condition(self, tag, field, val) :
def add_boundary_condition(self, tag, field, val = None, cb = None) :
if cb is not None :
val = 0.
vdof = self._dof_by_dimension[0]
edof = self._dof_by_dimension[self._fields[field][2].dimension()]
fixed_start = len(self._fixedidx)
fixedvid = []
for en in self._mesh.entities :
for ph in en.physicals :
phname = self._mesh.getPhysicalName(en.dimension, ph)
......@@ -84,12 +68,9 @@ class DofManager :
self._fixed[dof] = len(self._fixedval)
self._fixedval.append(val)
self._fixedidx.append(dof)
for dof in edof.get((field, e.tag), ()) :
if not dof in self._fixed :
self._fixed[dof] = len(self._fixedval)
self._fixedval.append(val)
self._fixedidx.append(dof)
fmap = np.array((self._ndof), np.int)
fixedvid.append(vid)
if cb is not None :
self._boundary_cb.append((cb, slice(fixed_start, len(self._fixedidx)), fixedvid))
class matrix(np.ndarray) :
......@@ -113,7 +94,7 @@ class DofManager :
class vector(np.ndarray) :
def __new__(cls, dofmanager) :
obj = np.zeros(dofmanager.getN()).view(cls)
obj = np.zeros(dofmanager.get_n()).view(cls)
obj._dofmanager = dofmanager
return obj
......@@ -132,7 +113,7 @@ class DofManager :
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(), 2))
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])
......@@ -150,7 +131,7 @@ class DofManager :
def gradient_at_points(self, els, xi, j) :
dofmanager = self._dofmanager
val = np.ndarray((len(els), dofmanager.n_fields(), 2))
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)
......@@ -159,19 +140,26 @@ class DofManager :
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)
NTot = sum((self.get_field_basis(i).n() for i in range(self.n_fields())))
self._jaccsrrows = np.ndarray((len(self.getElements(0)), NTot, NTot))
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]
......@@ -184,3 +172,4 @@ class DofManager :
self._fixedval = []
self._fixedidx = []
self._fixed = {}
self._boundary_cb = []
import numpy as np
class TriangleP1 :
points = np.array([[0., 0.], [1., 0.], [0., 1.]])
@staticmethod
def f(pts) :
r = np.ndarray((pts.shape[0], 3))
r[:, 0] = 1 - pts[:, 0] - pts[:, 1]
r[:, 1] = pts[:, 0]
r[:, 2] = pts[:, 1]
return r
@staticmethod
def df(pts) :
d = np.ndarray((pts.shape[0], 3, 2))
d[:, 0, 0] = -1
d[:, 0, 1] = -1
d[:, 1, 0] = 1
d[:, 1, 1] = 0
d[:, 2, 0] = 0
d[:, 2, 1] = 1
return d
@staticmethod
def dimension() :
return 2
@staticmethod
def n(dim=None) :
if dim is None :
return 3
elif dim == 0 :
return 1
else :
return 0
class TriangleP1Mini :
points = np.array([[0., 0.], [1., 0.], [0., 1.], [1./3, 1./3]])
@staticmethod
def f(pts) :
r = np.ndarray((pts.shape[0], 4))
......@@ -40,36 +79,93 @@ class TriangleP1Mini :
return 0
class TriangleP1 :
class TetrahedronP1 :
points = np.array([
[0., 0., 0.],
[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]
])
@staticmethod
def f(pts) :
r = np.ndarray((pts.shape[0], 3))
r[:, 0] = 1 - pts[:, 0] - pts[:, 1]
r = np.ndarray((pts.shape[0], 4))
r[:, 0] = 1 - pts[:, 0] - pts[:, 1] - pts[:, 2]
r[:, 1] = pts[:, 0]
r[:, 2] = pts[:, 1]
r[:, 3] = pts[:, 2]
return r
@staticmethod
def df(pts) :
d = np.ndarray((pts.shape[0], 3, 2))
d = np.zeros((pts.shape[0], 4, 3))
d[:, 0, 0] = -1
d[:, 0, 1] = -1
d[:, 0, 2] = -1
d[:, 1, 0] = 1
d[:, 1, 1] = 0
d[:, 2, 0] = 0
d[:, 2, 1] = 1
d[:, 3, 2] = 1
return d
@staticmethod
def dimension() :
return 2
return 3
@staticmethod
def n(dim=None) :
if dim is None :
return 3
return 4
elif dim == 0 :
return 1
else :
return 0
class TetrahedronP1Mini :
points = np.array([
[0., 0., 0.],
[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.],
[0.25, 0.25, 0.25]
])
@staticmethod
def f(pts) :
r = np.ndarray((pts.shape[0], 5))
r[:, 0] = 1 - pts[:, 0] - pts[:, 1] - pts[:, 2]
r[:, 1] = pts[:, 0]
r[:, 2] = pts[:, 1]
r[:, 3] = pts[:, 2]
r[:, 4] = r[:, 0] * r[:, 1] * r[:, 2] * r[:, 3]
return r
@staticmethod
def df(pts) :
d = np.zeros((pts.shape[0], 5, 3))
d[:, 0, 0] = -1
d[:, 0, 1] = -1
d[:, 0, 2] = -1
d[:, 1, 0] = 1
d[:, 2, 1] = 1
d[:, 3, 2] = 1
d[:, 4, 0] = pts[:, 1] * pts[:, 2] * (1 - 2*pts[:, 0] - pts[:, 1] - pts[:, 2])
d[:, 4, 1] = pts[:, 0] * pts[:, 2] * (1 - pts[:, 0] - 2*pts[:, 1] - pts[:, 2])
d[:, 4, 2] = pts[:, 0] * pts[:, 1] * (1 - pts[:, 0] - pts[:, 1] - 2*pts[:, 2])
return d
@staticmethod
def dimension() :
return 3
@staticmethod
def n(dim=None) :
if dim is None :
return 5
elif dim == 0 :
return 1
elif dim == 3 :
return 1
else :
return 0
......@@ -3,8 +3,6 @@ from . import legendre
import numpy as np
import scipy.sparse
D = 2
class MeshJacobian :
def update(self, x) :
self.x = x
......@@ -57,7 +55,7 @@ class MeshJacobian :
jfront = (els,) + (np.newaxis,) * (axis - 1) + (slice(None),)
jback = (np.newaxis,) * (naxes - axis - 1)
X = x[xfront + (0,) + xback] * self.dxidx[jfront + (0,) + jback]
for i in range(1,D) :
for i in range(1,self.dim) :
X += x[xfront + (i,) + xback] * self.dxidx[jfront + (i,) + jback]
return X
......@@ -68,7 +66,7 @@ class MeshJacobian :
jfront = (els,) + (np.newaxis,) * (axis - 1)
jback = (slice(None),) + (np.newaxis,) * (naxes - axis - 1)
X = x[xfront + (0,) + xback] * self.dxidx[jfront + (0,) + jback]
for i in range(1,D) :
for i in range(1,self.dim) :
X += x[xfront + (i,) + xback] * self.dxidx[jfront + (i,) + jback]
return X
......@@ -76,9 +74,12 @@ class MeshJacobian :
def __init__(self, mesh, dim) :
d = dofManager.DofManager(mesh)
self.dofManager = d
d.add_field(legendre.TriangleP1)
d.complete()
self.dim = dim
if self.dim == 3 :
d.add_field(legendre.TetrahedronP1)
else :
d.add_field(legendre.TriangleP1)
d.complete()
self.mesh = mesh
self.els = d.getElements(0)
self.update(np.array([[i[:dim] for i in el.vertices] for el in self.els]))
This diff is collapsed.
......@@ -352,6 +352,18 @@ static void _bboxtol(double *bbmin, double *bbmax) {
}
}
double volTet(double x0[3], double x1[3], double x2[3], double x3[3]) {
// x1 * (y2 * (z3 - z4) - z2 * (y3 - y4) + y3 * z4 - y4 * z3)
//- y1 * (x2 * (z3 - z4) - z2 * (x3 - x4) + x3 * z4 - x4 * z3)
//+ z1 * (x2 * (y3 - y4) - y2 * (x3 - x4) + x3 * y4 - x4 * y3)
//- (x2 * (y3 * z4 - y4 * z3) - y2 * (x3 * z4 - x4 * z3) + z2 * (x3 * y4 - x4 * y3))
return x0[0] * (x1[1] * (x2[2] - x3[2]) - x1[2] * (x2[1] - x3[1]) + x2[1] * x3[2] - x3[1] * x2[2])
- x0[1] * (x1[0] * (x2[2] - x3[2]) - x1[2] * (x2[0] - x3[0]) + x2[0] * x3[2] - x3[0] * x2[2])
+ x0[2] * (x1[0] * (x2[1] - x3[1]) - x1[1] * (x2[0] - x3[0]) + x2[0] * x3[1] - x3[0] * x2[1])
- (x1[0] * (x2[1] * x3[2] - x3[1] * x2[2]) - x1[1] * (x2[0] * x3[2] - x3[0] * x2[2]) + x1[2] * (x2[0] * x3[1] - x3[0] * x2[1]));
}
void ParticleToMesh(size_t nParticles, double *particles, size_t nElements, double *elements, size_t *elid, double *Xi)
{
double bbmin[DIMENSION], bbmax[DIMENSION];
......@@ -393,21 +405,44 @@ void ParticleToMesh(size_t nParticles, double *particles, size_t nElements, doub
cellSearch(tree, x, x, &found);
for (size_t j = 0; j < vectorSize(found); ++j) {
double *el = &elements[found[j] * N * DIMENSION];
double *X[3] = {el, el + DIMENSION, el + DIMENSION * 2};
double dx = x[0] - X[0][0], dy = x[1] - X[0][1];
double DX[2] = {X[1][0] - X[0][0], X[2][0] - X[0][0]};
double DY[2] = {X[1][1] - X[0][1], X[2][1] - X[0][1]};
double det = DX[1] * DY[0] - DX[0] * DY[1];
double xi = (DX[1] * dy - DY[1] * dx) / det;
double eta = -(DX[0] * dy - DY[0] * dx) / det;
if (xi > -1e-8 && eta > -1e-8 && 1 - xi - eta > -1e-8) {
Xi[i * DIMENSION + 0] = xi;
Xi[i * DIMENSION + 1] = eta;
if (DIMENSION == 2) {
double *X[3] = {el, el + DIMENSION, el + DIMENSION * 2};
double dx = x[0] - X[0][0], dy = x[1] - X[0][1];
double DX[2] = {X[1][0] - X[0][0], X[2][0] - X[0][0]};
double DY[2] = {X[1][1] - X[0][1], X[2][1] - X[0][1]};
double det = DX[1] * DY[0] - DX[0] * DY[1];
double xi = (DX[1] * dy - DY[1] * dx) / det;
double eta = -(DX[0] * dy - DY[0] * dx) / det;
if (xi > -1e-8 && eta > -1e-8 && 1 - xi - eta > -1e-8) {
Xi[i * DIMENSION + 0] = xi;
Xi[i * DIMENSION + 1] = eta;
elid[i] = found[j];
break;
}
}
else {
double *X[4] = {el, el + DIMENSION, el + DIMENSION * 2, el + DIMENSION * 3};
double v = volTet(X[0], X[1], X[2], X[3]);
double v0 = volTet(x, X[1], X[2], X[3]);
if (v0 * v < -1e-8)
continue;
double v1 = volTet(X[0], x, X[2], X[3]);
if (v1 * v < -1e-8)
continue;
double v2 = volTet(X[0], X[1], x, X[3]);
if (v2 * v < -1e-8)
continue;
double v3 = volTet(X[0], X[1], X[2], x);
if (v3 * v < -1e-8)
continue;
Xi[i * 3 + 0] = v1 / v;
Xi[i * 3 + 1] = v2 / v;
Xi[i * 3 + 2] = v3 / v;
elid[i] = found[j];
break;
}
}
if (elid[i] == -1) {
/*if (elid[i] == -1) {
printf(" PARTICLE %i OUTSIDE DOMAIN N2 search!!!\n", i);
double toll = -1000;
for (size_t j = 0; j < nElements; ++j) {
......@@ -426,12 +461,12 @@ void ParticleToMesh(size_t nParticles, double *particles, size_t nElements, doub
elid[i] = j;
break;
}
}
}*/
if (elid[i] == -1) {
printf(" PARTICLE %i OUTSIDE DOMAIN!!! %g %g (%g)\n", i, x[0], x[1], toll);
printf(" PARTICLE %i OUTSIDE DOMAIN!!! %g %g\n", i, x[0], x[1]);
exit(1);
}
}
//}
}
cellFree(tree);
}
......
#!/usr/bin/env python
import scontact2
import numpy as np
import mesh
import os
meshBoundary = True
outputdir = "output_depot"