Commit 93e7ae71 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

use material/boundary string names

parent d6ce95e8
Pipeline #4897 passed with stage
in 27 seconds
...@@ -34,6 +34,7 @@ def _write_array(f,anc,attr) : ...@@ -34,6 +34,7 @@ def _write_array(f,anc,attr) :
np.dtype('<i8'):b"Int64", np.dtype('<i8'):b"Int64",
np.dtype('<f4'):b"Float32", np.dtype('<f4'):b"Float32",
np.dtype('<f8'):b"Float64", np.dtype('<f8'):b"Float64",
np.dtype('uint8'):b"UInt8",
np.dtype('uint64'):b"UInt64" np.dtype('uint64'):b"UInt64"
}[a.dtype] }[a.dtype]
f.write(b'<DataArray %s type="%s" format="binary">\n'%(attr,tname)) f.write(b'<DataArray %s type="%s" format="binary">\n'%(attr,tname))
...@@ -44,7 +45,7 @@ def _write_array(f,anc,attr) : ...@@ -44,7 +45,7 @@ def _write_array(f,anc,attr) :
def _read_array(a) : def _read_array(a) :
raw = a.text.strip() raw = a.text.strip()
typename = a.attrib["type"] typename = a.attrib["type"]
dtype = {"Float64":np.float64,"Float32":np.float32,"Int64":np.int64,"Int32":np.int32,"UInt64":np.uint64}[typename] dtype = {"UInt8":np.uint8,"Float64":np.float64,"Float32":np.float32,"Int64":np.int64,"Int32":np.int32,"UInt64":np.uint64}[typename]
_,n,_,s = struct.unpack("iiii",b64decode(raw[:24])) _,n,_,s = struct.unpack("iiii",b64decode(raw[:24]))
data = zlib.decompress(b64decode(raw[24:24-(-s//3)*4])) data = zlib.decompress(b64decode(raw[24:24-(-s//3)*4]))
nc,nt = -1,-1 nc,nt = -1,-1
...@@ -196,3 +197,15 @@ def write_multipart(basename,parts,i,t,reset_index=None): ...@@ -196,3 +197,15 @@ def write_multipart(basename,parts,i,t,reset_index=None):
f.write(b'</vtkMultiBlockDataSet>\n'); f.write(b'</vtkMultiBlockDataSet>\n');
f.write(b'</VTKFile>\n'); f.write(b'</VTKFile>\n');
_write_index(basename+".pvd",os.path.basename(filename),i,t,reset_index) _write_index(basename+".pvd",os.path.basename(filename),i,t,reset_index)
def string_array_encode(a) :
joined = b"".join(a)
fmt = b"i%ii%is"%(len(a),len(joined))
r = struct.pack(fmt,len(a),*(len(i) for i in a),joined)
return np.array(list((i for i in r)),dtype="u1").reshape([-1,1])
def string_array_decode(e) :
n = struct.unpack("i",e[:4])[0]
l = struct.unpack("%ii"%n,e[4:4+4*n])
return struct.unpack("".join("%is"%i for i in l), e[4+n*4:])
...@@ -133,12 +133,19 @@ class ParticleProblem : ...@@ -133,12 +133,19 @@ class ParticleProblem :
def disk_tag(self): def disk_tag(self):
return self._get_idx("DiskTag") return self._get_idx("DiskTag")
def disk_material(self):
return self._get_idx("DiskMaterial")
def segments(self): def segments(self):
return self._get_matrix("Segment",3*self._dim+(1 if self._friction_enabled else 0)) return self._get_matrix("Segment",3*self._dim+(1 if self._friction_enabled else 0))
def segment_tag(self): def segment_tag(self):
return self._get_idx("SegmentTag") return self._get_idx("SegmentTag")
def segment_material(self):
return self._get_idx("SegmentMaterial")
def triangles(self): def triangles(self):
if self._dim == 3 : if self._dim == 3 :
return self._get_matrix("Triangle",12) return self._get_matrix("Triangle",12)
...@@ -151,6 +158,12 @@ class ParticleProblem : ...@@ -151,6 +158,12 @@ class ParticleProblem :
else : else :
return np.array([],np.int32) return np.array([],np.int32)
def triangle_material(self):
if self._dim == 3 :
return self._get_idx("TriangleMaterial")
else :
return np.array([],np.int32)
def set_boundary_velocity(self, tag, v) : def set_boundary_velocity(self, tag, v) :
self._lib.particleProblemGetTagId.restype = c_size_t self._lib.particleProblemGetTagId.restype = c_size_t
tag = self._lib.particleProblemGetTagId(self._p, tag.encode()) tag = self._lib.particleProblemGetTagId(self._p, tag.encode())
...@@ -173,10 +186,9 @@ class ParticleProblem : ...@@ -173,10 +186,9 @@ class ParticleProblem :
assert self._friction_enabled assert self._friction_enabled
self._lib.particleProblemSetFrictionRelaxation(self._p,c_double(relaxation)) self._lib.particleProblemSetFrictionRelaxation(self._p,c_double(relaxation))
def add_friction(self, static) : def set_friction_coefficient(self, mu, mat0="default",mat1="default") :
assert self._friction_enabled assert self._friction_enabled
arg_type = c_double*len(static) self._lib.particleProblemSetFrictionCoefficient(self._p, c_double(mu), mat0.encode(),mat1.encode())
self._lib.particleProblemAddFriction(self._p, arg_type(*static),len(static))
def iterate(self, dt, forces,tol=1e-8) : def iterate(self, dt, forces,tol=1e-8) :
"""Compute iteratively the set of velocities that respect the non-interpenetration constrain. """Compute iteratively the set of velocities that respect the non-interpenetration constrain.
...@@ -198,12 +210,15 @@ class ParticleProblem : ...@@ -198,12 +210,15 @@ class ParticleProblem :
else : else :
omega = np.zeros((v.shape[0],1)) omega = np.zeros((v.shape[0],1))
x = self.position() x = self.position()
tag = self._get_idx("ParticleTag").reshape(-1,1) material = self._get_idx("ParticleMaterial").reshape(-1,1)
if self._dim == 2 : if self._dim == 2 :
v = np.insert(v,self._dim,0,axis=1) v = np.insert(v,self._dim,0,axis=1)
x = np.insert(x,self._dim,0,axis=1) x = np.insert(x,self._dim,0,axis=1)
data = [("Mass",self.mass()), ("Radius",self.r()),("Velocity",v),("Omega",omega),("Tag",tag)] data = [("Mass",self.mass()), ("Radius",self.r()),("Velocity",v),("Omega",omega),("Material",material)]
VTK.write(odir+"/particles",i,t,None,x,data,vtktype="vtp") nmat = self._lib.particleProblemNMaterial(self._p)
self._lib.particleProblemGetMaterialTagName.restype = c_char_p
tags = list([self._lib.particleProblemGetMaterialTagName(self._p,i) for i in range(nmat)])
VTK.write(odir+"/particles",i,t,None,x,data,vtktype="vtp",field_data=[(b"MaterialNames",VTK.string_array_encode(tags))])
xdisk = self._get_matrix("Disk",2*self._dim+(2 if self._friction_enabled else 1))[:,:self._dim] xdisk = self._get_matrix("Disk",2*self._dim+(2 if self._friction_enabled else 1))[:,:self._dim]
nd = xdisk.shape[0] nd = xdisk.shape[0]
...@@ -217,8 +232,12 @@ class ParticleProblem : ...@@ -217,8 +232,12 @@ class ParticleProblem :
connectivity = np.arange(0,x.shape[0],dtype=np.int32) connectivity = np.arange(0,x.shape[0],dtype=np.int32)
types = np.repeat(np.array([1,3,5],dtype=np.int32),[nd,ns,nt]) types = np.repeat(np.array([1,3,5],dtype=np.int32),[nd,ns,nt])
offsets = np.cumsum(np.repeat([1,2,3],[nd,ns,nt]),dtype=np.int32) offsets = np.cumsum(np.repeat([1,2,3],[nd,ns,nt]),dtype=np.int32)
tags = np.array(np.hstack([self.disk_tag(),self.segment_tag(),self.triangle_tag()]),dtype=np.float64) tags = np.array(np.hstack([self.disk_tag(),self.segment_tag(),self.triangle_tag()]))
VTK.write(odir+"/boundaries",i,t,(types,connectivity,offsets),x,cell_data=[("Tag",tags.reshape([-1,1]))]) materials = np.array(np.hstack([self.disk_material(),self.segment_material(),self.triangle_material()]))
ntagname = self._lib.particleProblemNTag(self._p)
self._lib.particleProblemGetTagName.restype = c_char_p
tagnames = list([self._lib.particleProblemGetTagName(self._p,i) for i in range(ntagname)])
VTK.write(odir+"/boundaries",i,t,(types,connectivity,offsets),x,cell_data=[("Tag",tags.reshape([-1,1])),("Material",materials.reshape([-1,1]))],field_data=[(b"TagNames",VTK.string_array_encode(tagnames))])
self._lib.particleProblemContactN.restype = c_size_t self._lib.particleProblemContactN.restype = c_size_t
fdata = [] fdata = []
...@@ -238,31 +257,38 @@ class ParticleProblem : ...@@ -238,31 +257,38 @@ class ParticleProblem :
def read_vtk(self,dirname,i) : def read_vtk(self,dirname,i) :
"""Read an existing output file to set the particle data.""" """Read an existing output file to set the particle data."""
x,_,d,cdata,fdata = VTK.read(dirname+"/particles_%05d.vtp"%i) x,_,d,cdata,fdata = VTK.read(dirname+"/particles_%05d.vtp"%i)
matnames = VTK.string_array_decode(fdata["MaterialNames"])
matmap = {}
for j,n in enumerate(matnames) :
matmap[j] = self._lib.particleProblemGetMaterialTagId(self._p,n)
partmat =np.vectorize(matmap.get)(d["Material"])
self._lib.particleProblemClearParticles(self._p) self._lib.particleProblemClearParticles(self._p)
self._lib.particleProblemAddParticles(self._p,c_int(x.shape[0]), self._lib.particleProblemAddParticles(self._p,c_int(x.shape[0]),
_np2c(x[:,:self._dim]), _np2c(d["Mass"]), _np2c(d["Radius"]), _np2c(d["Velocity"][:,:self._dim]),_np2c(d["Omega"]),_np2c(d["Tag"])) _np2c(x[:,:self._dim]), _np2c(d["Mass"]), _np2c(d["Radius"]), _np2c(d["Velocity"][:,:self._dim]),_np2c(d["Omega"]),_np2c(partmat,dtype=np.int32))
x,cells,d,cdata,fdata = VTK.read(dirname+"/boundaries_%05d.vtu"%i) x,cells,d,cdata,fdata = VTK.read(dirname+"/boundaries_%05d.vtu"%i)
tagnames = VTK.string_array_decode(fdata["TagNames"])
tagmap = {}
for j,n in enumerate(tagnames) :
tagmap[j] = self._lib.particleProblemGetTagId(self._p,n)
self._lib.particleProblemClearBoundaries(self._p) self._lib.particleProblemClearBoundaries(self._p)
offsets = np.hstack([[0],cells["offsets"]]) offsets = np.hstack([[0],cells["offsets"]])
el = cells["connectivity"] el = cells["connectivity"]
tags = cdata["Tag"] tags = np.vectorize(tagmap.get)(cdata["Tag"])
materials = cdata["Material"]
for idx in np.where(cells["types"] == 1)[0] : for idx in np.where(cells["types"] == 1)[0] :
x0 = x[el[offsets[idx]],:self._dim] x0 = x[el[offsets[idx]],:self._dim]
t = int(tags[idx,0]) self._lib.particleProblemAddBoundaryDiskTagId(self._p, self._coord_type(*x0), c_double(0.), c_int(tags[idx,0]), c_int(materials[idx,0]))
self._lib.particleProblemAddBoundaryDiskTagId(self._p, self._coord_type(*x0), c_double(0.), c_int(t))
for idx in np.where(cells["types"] == 3)[0] : for idx in np.where(cells["types"] == 3)[0] :
x0 = x[el[offsets[idx]],:self._dim] x0 = x[el[offsets[idx]],:self._dim]
x1 = x[el[offsets[idx]+1],:self._dim] x1 = x[el[offsets[idx]+1],:self._dim]
t = int(tags[idx,0]) self._lib.particleProblemAddBoundarySegmentTagId(self._p, self._coord_type(*x0), self._coord_type(*x1), c_int(tags[idx,0]), c_int(materials[idx,0]))
self._lib.particleProblemAddBoundarySegmentTagId(self._p, self._coord_type(*x0), self._coord_type(*x1), c_int(t))
for idx in np.where(cells["types"] == 5)[0] : for idx in np.where(cells["types"] == 5)[0] :
x0 = x[el[offsets[idx]],:self._dim] x0 = x[el[offsets[idx]],:self._dim]
x1 = x[el[offsets[idx]+1],:self._dim] x1 = x[el[offsets[idx]+1],:self._dim]
x2 = x[el[offsets[idx]+2],:self._dim] x2 = x[el[offsets[idx]+2],:self._dim]
t = int(tags[idx,0]) self._lib.particleProblemAddBoundaryTriangleTagId(self._p, self._coord_type(*x0), self._coord_type(*x1), self._coord_type(*x2), c_int(tags[idx,0]), c_int(materials[idx,0]))
self._lib.particleProblemAddBoundaryTriangleTagId(self._p, self._coord_type(*x0), self._coord_type(*x1), self._coord_type(*x2), c_int(t))
_,_,_,_,fdata = VTK.read(dirname+"/contacts_%05d.vtp"%i) _,_,_,_,fdata = VTK.read(dirname+"/contacts_%05d.vtp"%i)
ks = ("particle_particle","particle_disk","particle_segment","particle_triangle") ks = ("particle_particle","particle_disk","particle_segment","particle_triangle")
v = np.vstack([fdata[k] for k in ks]) v = np.vstack([fdata[k] for k in ks])
...@@ -272,19 +298,19 @@ class ParticleProblem : ...@@ -272,19 +298,19 @@ class ParticleProblem :
t = np.repeat([0,1,2,3],[fdata[k].shape[0] for k in ks]) t = np.repeat([0,1,2,3],[fdata[k].shape[0] for k in ks])
self._lib.particleProblemSetContacts(self._p,c_int(t.shape[0]),_np2c(o,np.uint64),_np2c(v),_np2c(t,np.int32)) self._lib.particleProblemSetContacts(self._p,c_int(t.shape[0]),_np2c(o,np.uint64),_np2c(v),_np2c(t,np.int32))
def add_boundary_disk(self, x0, r, tag) : def add_boundary_disk(self, x0, r, tag, material="default") :
self._lib.particleProblemAddBoundaryDisk.restype = c_size_t self._lib.particleProblemAddBoundaryDisk.restype = c_size_t
return self._lib.particleProblemAddBoundaryDisk(self._p, self._coord_type(*x0), c_double(r), tag.encode()) return self._lib.particleProblemAddBoundaryDisk(self._p, self._coord_type(*x0), c_double(r), tag.encode(),material.encode())
def add_boundary_segment(self, x0, x1, tag) : def add_boundary_segment(self, x0, x1, tag, material="default") :
self._lib.particleProblemAddBoundarySegment.restype = c_size_t self._lib.particleProblemAddBoundarySegment.restype = c_size_t
return self._lib.particleProblemAddBoundarySegment(self._p, self._coord_type(*x0), self._coord_type(*x1), tag.encode()) return self._lib.particleProblemAddBoundarySegment(self._p, self._coord_type(*x0), self._coord_type(*x1), tag.encode(),material.encode())
def add_boundary_triangle(self, x0, x1, x2, tag) : def add_boundary_triangle(self, x0, x1, x2, tag, material="default") :
if self._dim != 3 : if self._dim != 3 :
raise NameError("Triangle boundaries only available in 3D") raise NameError("Triangle boundaries only available in 3D")
self._lib.particleProblemAddBoundaryTriangle.restype = c_size_t self._lib.particleProblemAddBoundaryTriangle.restype = c_size_t
return self._lib.particleProblemAddBoundaryTriangle(self._p, self._coord_type(*x0), self._coord_type(*x1), self._coord_type(*x2), tag.encode()) return self._lib.particleProblemAddBoundaryTriangle(self._p, self._coord_type(*x0), self._coord_type(*x1), self._coord_type(*x2), tag.encode(),material.encode())
def add_particle(self, x, r, m, tag): def add_particle(self, x, r, m, tag):
"""Add a particle in the particle container. """Add a particle in the particle container.
...@@ -295,9 +321,9 @@ class ParticleProblem : ...@@ -295,9 +321,9 @@ class ParticleProblem :
m -- particle mass m -- particle mass
tag -- particle tag tag -- particle tag
""" """
self._lib.particleProblemAddParticle(self._p, self._coord_type(*x), c_double(r), c_double(m), c_int(tag)) self._lib.particleProblemAddParticle(self._p, self._coord_type(*x), c_double(r), c_double(m), tag.encode())
def load_msh_boundaries(self, filename, tags, shift=None) : def load_msh_boundaries(self, filename, tags, shift=None,material="default") :
"""Load the numerical domain and set the physical boundaries the particles cannot go through. """Load the numerical domain and set the physical boundaries the particles cannot go through.
Keyword arguments: Keyword arguments:
...@@ -321,10 +347,10 @@ class ParticleProblem : ...@@ -321,10 +347,10 @@ class ParticleProblem :
if v[3] in self._addv : if v[3] in self._addv :
continue continue
self._addv.add(v[3]) self._addv.add(v[3])
self.add_boundary_disk((v[0] + shift[0], v[1] + shift[1]), 0.000, stag) self.add_boundary_disk((v[0] + shift[0], v[1] + shift[1]), 0.000, stag,material)
v0 = el.vertices[1] v0 = el.vertices[1]
v1 = el.vertices[0] v1 = el.vertices[0]
self.add_boundary_segment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), stag) self.add_boundary_segment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), stag,material)
else : else :
self._adde = set() self._adde = set()
pnum = [m.getPhysicalNumber(2, i) for i in tags] pnum = [m.getPhysicalNumber(2, i) for i in tags]
...@@ -338,14 +364,14 @@ class ParticleProblem : ...@@ -338,14 +364,14 @@ class ParticleProblem :
v0 = el.vertices[j] v0 = el.vertices[j]
if not v0[3] in self._addv : if not v0[3] in self._addv :
self._addv.add(v0[3]) self._addv.add(v0[3])
self.add_boundary_disk((v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]),0.,stag) self.add_boundary_disk((v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]),0.,stag,material)
v1 = el.vertices[(j+1)%3] v1 = el.vertices[(j+1)%3]
if (not (v0[3],v1[3]) in self._adde) and (not (v1[3],v0[3]) in self._adde) : if (not (v0[3],v1[3]) in self._adde) and (not (v1[3],v0[3]) in self._adde) :
self._adde.add((v0[3],v1[3])) self._adde.add((v0[3],v1[3]))
self.add_boundary_segment( self.add_boundary_segment(
(v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]), (v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]),
(v1[0]+shift[0],v1[1]+shift[1],v1[2]+shift[2]), (v1[0]+shift[0],v1[1]+shift[1],v1[2]+shift[2]),
stag) stag,material)
v0 = el.vertices[1] v0 = el.vertices[1]
v1 = el.vertices[0] v1 = el.vertices[0]
...@@ -354,4 +380,4 @@ class ParticleProblem : ...@@ -354,4 +380,4 @@ class ParticleProblem :
(v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]), (v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]),
(v1[0]+shift[0],v1[1]+shift[1],v1[2]+shift[2]), (v1[0]+shift[0],v1[1]+shift[1],v1[2]+shift[2]),
(v2[0]+shift[0],v2[1]+shift[1],v2[2]+shift[2]), (v2[0]+shift[0],v2[1]+shift[1],v2[2]+shift[2]),
stag) stag,material)
...@@ -42,8 +42,10 @@ struct _ParticleProblem{ ...@@ -42,8 +42,10 @@ struct _ParticleProblem{
Contact *contacts; Contact *contacts;
double *position, *velocity, *externalForces; double *position, *velocity, *externalForces;
Disk *disks; Disk *disks;
char **tagname; char **tagname, **materialName;
int *diskTag, *segmentTag, *particleTag; int *diskTag, *segmentTag;
int *diskMaterial, *segmentMaterial;
int *particleMaterial;
Segment *segments; Segment *segments;
#if FRICTION_ENABLED #if FRICTION_ENABLED
double *omega; double *omega;
...@@ -53,6 +55,7 @@ struct _ParticleProblem{ ...@@ -53,6 +55,7 @@ struct _ParticleProblem{
#if DIMENSION == 3 #if DIMENSION == 3
Triangle *triangles; Triangle *triangles;
int *triangleTag; int *triangleTag;
int *triangleMaterial;
#endif #endif
int use_queue; int use_queue;
}; };
...@@ -138,9 +141,11 @@ static void diskBoundingBox(const Disk *d, double *pmin, double *pmax) { ...@@ -138,9 +141,11 @@ static void diskBoundingBox(const Disk *d, double *pmin, double *pmax) {
} }
} }
size_t particleProblemAddBoundaryDiskTagId(ParticleProblem *p, const double x[DIMENSION], double r, size_t tag) {
size_t particleProblemAddBoundaryDiskTagId(ParticleProblem *p, const double x[DIMENSION], double r, int tag, int materialTag) {
Disk *d = vectorPush(&p->disks); Disk *d = vectorPush(&p->disks);
*vectorPush(&p->diskTag) = tag; *vectorPush(&p->diskTag) = tag;
*vectorPush(&p->diskMaterial) = materialTag;
d->r = r; d->r = r;
#ifdef FRICTION_ENABLED #ifdef FRICTION_ENABLED
d->vt = 0; d->vt = 0;
...@@ -152,8 +157,8 @@ size_t particleProblemAddBoundaryDiskTagId(ParticleProblem *p, const double x[DI ...@@ -152,8 +157,8 @@ size_t particleProblemAddBoundaryDiskTagId(ParticleProblem *p, const double x[DI
return vectorSize(p->disks) - 1; return vectorSize(p->disks) - 1;
} }
size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x[DIMENSION], double r, const char *tag) { size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x[DIMENSION], double r, const char *tag, const char *materialTag) {
return particleProblemAddBoundaryDiskTagId(p, x, r, particleProblemGetTagId(p, tag)); return particleProblemAddBoundaryDiskTagId(p,x,r,particleProblemGetTagId(p,tag),particleProblemGetMaterialTagId(p,materialTag));
} }
static int diskInitContact(size_t id, const Disk *d, size_t particle, const Particle *p, double *x, double alert, Contact *c) { static int diskInitContact(size_t id, const Disk *d, size_t particle, const Particle *p, double *x, double alert, Contact *c) {
...@@ -255,9 +260,10 @@ struct _Triangle { ...@@ -255,9 +260,10 @@ struct _Triangle {
double v[DIMENSION]; double v[DIMENSION];
}; };
void particleProblemAddBoundaryTriangleTagId(ParticleProblem *p, const double x0[3], const double x1[3], const double x2[3], size_t tag) { void particleProblemAddBoundaryTriangleTagId(ParticleProblem *p, const double x0[3], const double x1[3], const double x2[3], int tag, int materialTag) {
Triangle *t = vectorPush(&p->triangles); Triangle *t = vectorPush(&p->triangles);
*vectorPush(&p->triangleTag) = tag; *vectorPush(&p->triangleTag) = tag;
*vectorPush(&p->triangleMaterial) = materialTag;
for (int i = 0; i < DIMENSION; ++i) { for (int i = 0; i < DIMENSION; ++i) {
t->p[0][i] = x0[i]; t->p[0][i] = x0[i];
t->p[1][i] = x1[i]; t->p[1][i] = x1[i];
...@@ -266,8 +272,8 @@ void particleProblemAddBoundaryTriangleTagId(ParticleProblem *p, const double x0 ...@@ -266,8 +272,8 @@ void particleProblemAddBoundaryTriangleTagId(ParticleProblem *p, const double x0
} }
} }
void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[3], const double x1[3], const double x2[3], const char *tag) { void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[3], const double x1[3], const double x2[3], const char *tag, const char *material) {
return particleProblemAddBoundaryTriangleTagId(p, x0, x1, x2, particleProblemGetTagId(p, tag)); return particleProblemAddBoundaryTriangleTagId(p, x0, x1, x2, particleProblemGetTagId(p, tag), particleProblemGetMaterialTagId(p, material));
} }
static void triangleBoundingBox(const Triangle *t, double *pmin, double *pmax) { static void triangleBoundingBox(const Triangle *t, double *pmin, double *pmax) {
...@@ -340,6 +346,12 @@ static int triangleInitContact(size_t id, const Triangle *t, size_t particle, co ...@@ -340,6 +346,12 @@ static int triangleInitContact(size_t id, const Triangle *t, size_t particle, co
#endif #endif
#ifdef FRICTION_ENABLED
static double particleProblemGetMu(const ParticleProblem *p, int mat0, int mat1) {
return p->mu[mat0*particleProblemNMaterial(p)+mat1];
}
#endif
static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double dt, double tol) { static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double dt, double tol) {
double vn = c->dv; double vn = c->dv;
for (int i = 0; i<DIMENSION; ++i) for (int i = 0; i<DIMENSION; ++i)
...@@ -357,8 +369,7 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d ...@@ -357,8 +369,7 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
(p->velocity[c->o0*2+1]-p->velocity[c->o1*2+1])*t[1]+ (p->velocity[c->o0*2+1]-p->velocity[c->o1*2+1])*t[1]+
p->omega[c->o0]*p->particles[c->o0].r + p->omega[c->o0]*p->particles[c->o0].r +
p->omega[c->o1]*p->particles[c->o1].r; p->omega[c->o1]*p->particles[c->o1].r;
int index = p->particleTag[c->o0]*vectorSize(p->mu) + p->particleTag[c->o1]; const double mu = particleProblemGetMu(p,p->particleMaterial[c->o0],p->particleMaterial[c->o1]);
double mu = p->mu[index];
if(vt>(7./2.)*mu*(dp)){ if(vt>(7./2.)*mu*(dp)){
ct = p->friction_relaxation*mu*(dp); ct = p->friction_relaxation*mu*(dp);
} }
...@@ -410,8 +421,7 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double dt, d ...@@ -410,8 +421,7 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double dt, d
p->velocity[c->o1 * DIMENSION+1]*t[1]+ p->velocity[c->o1 * DIMENSION+1]*t[1]+
p->omega[c->o1]*p->particles[c->o1].r+ p->omega[c->o1]*p->particles[c->o1].r+
p->disks[c->o0].vt; p->disks[c->o0].vt;
int index = (p->particleTag[c->o1]+1)*vectorSize(p->mu) -1; const double mu = particleProblemGetMu(p,p->particleMaterial[c->o1],p->diskMaterial[c->o0]);
double mu = p->mu[index];
double ct = 0; double ct = 0;
if(vt>(7./2.)*mu*dp){ if(vt>(7./2.)*mu*dp){
ct = p->friction_relaxation*mu*dp; ct = p->friction_relaxation*mu*dp;
...@@ -457,8 +467,7 @@ static int contactParticleSegmentSolve(Contact *c, ParticleProblem *p, double dt ...@@ -457,8 +467,7 @@ static int contactParticleSegmentSolve(Contact *c, ParticleProblem *p, double dt
p->velocity[c->o1*DIMENSION +1]*t[1]+ p->velocity[c->o1*DIMENSION +1]*t[1]+
p->omega[c->o1]*p->particles[c->o1].r+ p->omega[c->o1]*p->particles[c->o1].r+
p->segments[c->o0].vt; p->segments[c->o0].vt;
int index = (p->particleTag[c->o1]+1)*vectorSize(p->mu) -1; const double mu = particleProblemGetMu(p,p->particleMaterial[c->o1],p->segmentMaterial[c->o0]);
double mu = p->mu[index];
double ct = 0; double ct = 0;
if(vt>(7./2.)*mu*dp){ if(vt>(7./2.)*mu*dp){
ct = p->friction_relaxation*mu*dp; ct = p->friction_relaxation*mu*dp;
...@@ -838,7 +847,6 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert) ...@@ -838,7 +847,6 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
t[0] = c->n[1]; t[0] = c->n[1];
t[1] = - c->n[0]; t[1] = - c->n[0];
c->ct = cold->ct; c->ct = cold->ct;
int index = p->particleTag[c->o0]*vectorSize(p->mu) + p->particleTag[c->o1];
coordAdd(&p->velocity[c->o0 * DIMENSION], -c->ct*c->a0, t); coordAdd(&p->velocity[c->o0 * DIMENSION], -c->ct*c->a0, t);
coordAdd(&p->velocity[c->o1 * DIMENSION], c->ct*c->a1, t); coordAdd(&p->velocity[c->o1 * DIMENSION], c->ct*c->a1, t);
p->omega[c->o0] -= c->In0*(c->ct); p->omega[c->o0] -= c->In0*(c->ct);
...@@ -867,7 +875,6 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert) ...@@ -867,7 +875,6 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
t[0] = c->n[1]; t[0] = c->n[1];
t[1] = - c->n[0]; t[1] = - c->n[0];
c->ct = cold->ct; c->ct = cold->ct;
int index = (p->particleTag[c->o1]+1)*vectorSize(p->mu) -1;
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*c->a1, t); coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*c->a1, t);
p->omega[c->o1] -= c->In1*(c->ct); p->omega[c->o1] -= c->In1*(c->ct);
#endif #endif
...@@ -892,7 +899,6 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert) ...@@ -892,7 +899,6 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
t[0] = c->n[1]; t[0] = c->n[1];
t[1] = -c->n[0]; t[1] = -c->n[0];
c->ct = cold->ct; c->ct = cold->ct;
int index = (p->particleTag[c->o1]+1)*vectorSize(p->mu)-1;
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*c->a1, t); coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*c->a1, t);
p->omega[c->o1] -= c->In1*(c->ct); p->omega[c->o1] -= c->In1*(c->ct);
#endif #endif
...@@ -1110,9 +1116,10 @@ void particleProblemIterate(ParticleProblem *p, double alert, double dt, double ...@@ -1110,9 +1116,10 @@ void particleProblemIterate(ParticleProblem *p, double alert, double dt, double
#endif #endif
} }
size_t particleProblemAddBoundarySegmentTagId(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], size_t tag) { size_t particleProblemAddBoundarySegmentTagId(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], int tag, int materialTag) {
Segment *s = vectorPush(&p->segments); Segment *s = vectorPush(&p->segments);
*vectorPush(&p->segmentTag) = tag; *vectorPush(&p->segmentTag) = tag;
*vectorPush(&p->segmentMaterial) = materialTag;
for (int i = 0; i < DIMENSION; ++i) { for (int i = 0; i < DIMENSION; ++i) {
s->p[0][i] = x0[i]; s->p[0][i] = x0[i];
s->p[1][i] = x1[i]; s->p[1][i] = x1[i];
...@@ -1124,11 +1131,11 @@ size_t particleProblemAddBoundarySegmentTagId(ParticleProblem *p, const double x ...@@ -1124,11 +1131,11 @@ size_t particleProblemAddBoundarySegmentTagId(ParticleProblem *p, const double x
return vectorSize(p->segments) - 1; return vectorSize(p->segments) - 1;
} }
size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const char *tag) { size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const char *tag, const char *material) {
return particleProblemAddBoundarySegmentTagId(p, x0, x1, particleProblemGetTagId(p, tag)); return particleProblemAddBoundarySegmentTagId(p, x0, x1, particleProblemGetTagId(p, tag),particleProblemGetMaterialTagId(p, material));
} }
void particleProblemAddParticle(ParticleProblem *p, const double x[DIMENSION], double r, double m, int tag) { void particleProblemAddParticleMaterialTagId(ParticleProblem *p, const double x[DIMENSION], double r, double m, int tag) {
size_t n = vectorSize(p->particles); size_t n = vectorSize(p->particles);
Particle *particle = vectorPush(&p->particles); Particle *particle = vectorPush(&p->particles