https://github.com/dbukenberger/HexahedralLattice
Tip revision: 991ca6893ac1feec28f491429111e0fbdcefccc0 authored by dbukenberger on 20 November 2024, 07:08:57 UTC
Update requirements.txt
Update requirements.txt
Tip revision: 991ca68
CubeObject.py
from femUtil import *
from PadObject import *
# quad and hex filter kernels
qKernels = [[1,-1,0,-1,1,0,0,0,0],
[0,-1,1,0,1,-1,0,0,0],
[0,0,0,-1,1,0,1,-1,0],
[0,0,0,0,1,-1,0,-1,1]]
hKernels = [[0,1,0,0,-1,0,0,0,0,0,-1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,1,-1,0,0,0,0,0,0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,-1,1,0,0,0,0,0,0,0,1,-1,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,-1,0,0,1,0,0,0,0,0,1,0,0,-1,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,1,-1,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,-1,1,0,1,-1,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,-1,1,0,1,-1,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,1,-1,0,-1,1,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,-1,0,0,1,0,0,0,0,0,1,0,0,-1,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,-1,1,0,0,0,0,0,0,0,1,-1,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,1,-1,0,0,0,0,0,0,0,-1,1,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,-1,0,0,0,0,0,-1,0,0,1,0],
[1,-1,0,-1,-1,0,0,0,0,-1,-1,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,-1,1,0,-1,-1,0,0,0,0,-1,-1,0,1,-1,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,-1,-1,0,1,-1,0,0,0,0,-1,1,0,-1,-1,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,-1,-1,0,-1,1,0,0,0,0,1,-1,0,-1,-1,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,-1,-1,0,-1,1,0,0,0,0,1,-1,0,-1,-1,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,-1,-1,0,1,-1,0,0,0,0,-1,1,0,-1,-1,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,-1,1,0,-1,-1,0,0,0,0,-1,-1,0,1,-1,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,1,-1,0,-1,-1,0,0,0,0,-1,-1,0,-1,1]]
mKernels = {2: qKernels, 3: hKernels}
class ConfigParams:
def __init__(self, objName, pLambda, lambdaRamp, autoScale, optiAlign, useNormals, stressWeighted, useConstraints, maxIters, stopDelta, gridRes, finalLayer, useProjection, useTransform, hullTolerance, hullSmoothing, autoSave = True):
self.objName = os.path.basename(objName).split('.')[0]
self.fileName = 'data/%s.stress'%self.objName
self.pLambda = pLambda
self.lambdaRamp = lambdaRamp
if self.lambdaRamp:
self.lambdaFun = lambda i: (i/self.lambdaRamp if i < self.lambdaRamp else 1) * self.pLambda
else:
self.lambdaFun = lambda i: self.pLambda
self.autoScale = autoScale
self.optiAlign = optiAlign
self.useNormals = useNormals
self.stressWeighted = stressWeighted
self.useConstraints = useConstraints
self.maxIters = maxIters
self.stopDelta = stopDelta
self.gridRes = gridRes
self.finalLayer = finalLayer
self.useProjection = useProjection
self.useTransform = useTransform
self.hullTolerance = hullTolerance
self.hullSmoothing = hullSmoothing
self.reset()
if autoSave:
self.save()
@classmethod
def fromFile(cls, fileName):
lines = np.loadtxt(fileName, str)
for line in lines:
if 'objName' in line[0]:
fileName = 'data/%s.stress'%line[1]
if 'pLambda' in line[0]:
pLambda = float(line[1])
if 'lambdaRamp' in line[0]:
lambdaRamp = int(line[1])
if 'autoScale' in line[0]:
autoScale = 'True' in line[1]
if 'optiAlign' in line[0]:
if 'False' in line[1]:
optiAlign = False
elif 'True' in line[1]:
optiAlign = True
else:
optiAlign = list(map(float, line[1].split(',')))
if 'useNormals' in line[0]:
useNormals = 'True' in line[1]
if 'stressWeighted' in line[0]:
stressWeighted = 'True' in line[1]
if 'useConstraints' in line[0]:
useConstraints = 'True' in line[1]
if 'maxIters' in line[0]:
maxIters = int(line[1])
if 'stopDelta' in line[0]:
stopDelta = float(line[1])
if 'gridRes' in line[0]:
gridRes = int(line[1])
if 'finalLayer' in line[0]:
finalLayer = 'True' in line[1]
if 'useProjection' in line[0]:
useProjection = 'True' in line[1]
if 'useTransform' in line[0]:
useTransform = 'True' in line[1]
if 'hullTolerance' in line[0]:
hullTolerance = 'True' in line[1]
if 'hullSmoothing' in line[0]:
if 'False' in line[1]:
hullSmoothing = False
elif 'True' in line[1]:
hullSmoothing = True
else:
hullSmoothing = float(line[1])
return cls(fileName, pLambda, lambdaRamp, autoScale, optiAlign, useNormals, stressWeighted, useConstraints, maxIters, stopDelta, gridRes, finalLayer, useProjection, useTransform, hullTolerance, hullSmoothing, False)
def reset(self):
self.rhoInit = 1e-3
self.abstol = 1e-5
self.reltol = 1e-3
self.mu = 10.0
self.tao = 2.0
self.maxIterADMM = 500
def save(self):
with open('data/%s.cfg'%self.objName, 'w') as fh:
fh.write('objName:\t\t%s\n'%self.objName)
fh.write('pLambda:\t\t%d\n'%self.pLambda)
fh.write('lambdaRamp:\t\t%d\n'%self.lambdaRamp)
fh.write('autoScale:\t\t%s\n'%self.autoScale)
if type(self.optiAlign) == bool:
fh.write('optiAlign:\t\t%s\n'%self.optiAlign)
else:
fh.write('optiAlign:\t\t'+(','.join(['%g']*len(self.optiAlign)))%tuple(self.optiAlign)+'\n')
fh.write('useNormals:\t\t%s\n'%self.useNormals)
fh.write('stressWeighted:\t%s\n'%self.stressWeighted)
fh.write('useConstraints:\t%s\n'%self.useConstraints)
fh.write('maxIters:\t\t%d\n'%self.maxIters)
fh.write('stopDelta:\t\t%g\n'%self.stopDelta)
fh.write('gridRes:\t\t%d\n'%self.gridRes)
fh.write('finalLayer:\t\t%s\n'%self.finalLayer)
fh.write('useProjection:\t%s\n'%self.useProjection)
fh.write('useTransform:\t%s\n'%self.useTransform)
fh.write('hullTolerance:\t%s\n'%self.hullTolerance)
fh.write('hullSmoothing:\t%s\n'%self.hullSmoothing)
class CubeObject:
def __init__(self, cfgParams):
self.cfgParams = cfgParams
if self.cfgParams.fileName.split('.')[-1] == 'stress':
self.verts, self.ts, self.flxIdxs, forceVecs, self.fixIdxs, self.vmStress, self.pStress, self.pStressE, self.sMats = loadStressFile(self.cfgParams.fileName)
self.aIdxs = np.concatenate([self.flxIdxs, self.fixIdxs])
self.nDim = self.verts.shape[1]
else:
return
if self.nDim == 2:
self.edges = facesToEdges(self.ts, False)
unq, inv, cnt = np.unique(cantorPiKV(self.edges), return_inverse = True, return_counts = True)
self.boundaryEdgeMask = cnt[inv] == 1
if self.cfgParams.useNormals:
self.boundaryVertMask = np.zeros(len(self.verts), np.bool_)
self.boundaryVertMask[np.unique(self.edges[self.boundaryEdgeMask].ravel())] = True
boundaryEdges = self.edges[self.boundaryEdgeMask]
boundaryLoop = edgesToPath(boundaryEdges)
bEdges = np.transpose([boundaryLoop, np.roll(boundaryLoop, 1)])
bEdgeVecs = self.verts[bEdges[:,1]] - self.verts[bEdges[:,0]]
bVertDirs = np.dot(bEdgeVecs, Mr2D(-np.pi/2)) + np.dot(np.roll(bEdgeVecs, -1, axis=0), Mr2D(-np.pi/2))
self.vNormals = np.zeros_like(self.verts)
self.vNormals[boundaryLoop] = normVec(bVertDirs)
self.pStress[self.aIdxs] = np.dstack([self.vNormals[self.aIdxs], np.dot(self.vNormals[self.aIdxs], Mr2D(np.pi/2))])
else:
self.tris = tetrasToFaces(self.ts)
unq, inv, cnt = np.unique(cantorPiKV(self.tris), return_inverse = True, return_counts = True)
self.boundaryTrisMask = cnt[inv] == 1
if self.cfgParams.useNormals:
if True: # default, normals only
self.vNormals = igl.per_vertex_normals(self.verts, self.tris[self.boundaryTrisMask])
m = np.bitwise_not(np.any(np.isnan(self.vNormals),axis=1))
aMask = np.zeros_like(m)
aMask[self.aIdxs] = True
for aIdx in self.aIdxs:
self.pStress[aIdx] = rotateAsToB(self.pStress[aIdx], self.vNormals[aIdx], self.pStress[aIdx,0])
else: # experimental, curvature tensor
ts = self.tris[self.boundaryTrisMask]
boundaryVertIdxs = np.unique(ts.ravel())
pc = igl.principal_curvature(self.verts[boundaryVertIdxs], reIndexIndices(ts))
allVertPc = np.zeros_like(self.pStress)
allVertPc[boundaryVertIdxs] = np.transpose(np.dstack([cross(pc[0], pc[1]), pc[0], pc[1]]), axes=(0,2,1))
self.pStress[self.aIdxs] = allVertPc[self.aIdxs]
print("#v % 6d, #e % 6d"%(len(self.verts), len(self.ts)))
self.isRestored = True
# scale object
self.BBinit = np.float32([self.verts.min(axis=0), self.verts.max(axis=0)])
if self.cfgParams.autoScale:
self.verts -= self.BBinit[0]
self.verts /= max(self.BBinit[1]-self.BBinit[0])/2
self.verts -= self.verts.max(axis=0)/2
self.isRestored = False
# orient object
self.MinitRot = np.eye(self.nDim)
if self.cfgParams.optiAlign:
if type(self.cfgParams.optiAlign) == list:
self.MinitRot = Mr2D(np.pi * self.cfgParams.optiAlign[0]) if self.nDim == 2 else Mr3D(np.pi * self.cfgParams.optiAlign[0], np.pi * self.cfgParams.optiAlign[1], np.pi * self.cfgParams.optiAlign[2])
else:
w = self.pStressE.sum(axis=1).reshape(-1,1,1) #* 0 + 1
self.MinitRot = computePrincipalStress((self.sMats * w/w.sum()).sum(axis=0))[0]
self.MinitRot = computeMinTransformation(self.MinitRot)
self.verts = np.dot(self.verts, self.MinitRot.T)
self.pStress = innerNxM(self.pStress, [self.MinitRot])
self.isRestored = False
self.sf = norm(self.verts.max(axis=0)-self.verts.min(axis=0)) / 1000
self.cVerts = self.verts.copy()
self.cVertss = [self.cVerts.copy()]
self.ks = []
self.energy = 0.0
self.relativeDeltaV = sys.float_info.max
self.tLocal = []
self.tGlobal = []
self.setupMatrices()
self.Rs = [self.R.copy()]
if self.cfgParams.useConstraints:
cIdxs = self.aIdxs
self.cRHS = np.zeros_like(self.verts)
self.cRHS[cIdxs] = self.verts[cIdxs]
E = sparse.dok_matrix((self.nVerts, self.nVerts), dtype=np.float64)
E[cIdxs,cIdxs] = 1
E = E.tocsc()
M = (self.L.T @ self.L + E).tocsc()
else:
M = (self.L.T @ self.L).tocsc()
# check sparseCond(M) and add eps * I
c = 1e-12
trys = 0
print('Construct solver', M.shape)
while not hasattr(self, 'solver'):
try:
Mplus = (M + sparse.eye(M.shape[0]).tocsc() * c) if trys else M
self.solver = pypardiso.factorized(Mplus) if pypardisoFound else spla.factorized(Mplus)
except RuntimeError as e:
if 'singular' in e.args[0]:
c *= 10
trys += 1
print('Desingularize trys:', trys)
if trys:
print('Mat singular fix: ', c, trys)
else:
print('SPLA solver done.')
def setupMatrices(self):
self.L = igl.cotmatrix(np.float64(self.verts), self.ts)
self.M = igl.massmatrix(self.verts, self.ts)
self.vMass = self.M.diagonal()
self.volumeInit = self.vMass.sum()
self.centerInit = (self.verts.min(axis=0) + self.verts.max(axis=0))/2
self.nVerts = len(self.verts)
self.nNeighbors = self.L.getnnz(axis=0).max()-1
nFan = 1 if self.nDim == 3 else 3
self.nIdxs = np.ones((self.nVerts,2,self.nNeighbors*nFan),np.int32) * np.arange(self.nVerts)[...,None,None]
self.nWeights = np.zeros_like(self.nIdxs[:,0], dtype=np.float32)
self.Ralpha = np.ones((self.nVerts, self.nNeighbors*nFan)) / 2
if self.nDim == 2:
tm = trimesh.Trimesh(self.verts, self.ts)
bLoop = np.bool_(igl.is_border_vertex(self.verts, self.ts))
for idx, col in tqdm(enumerate(self.L), total = self.L.shape[0], ascii=True, desc='weights'):
msk = col.indices != idx
indices = col.indices[msk]
values = col.data[msk]
if True or self.nDim == 3: # default, only spokes
self.nIdxs[idx,0,:len(indices)] = indices
self.nWeights[idx,:len(indices)] = values
else: # experimental for 2D, with rims
fIdxs = tm.vertex_faces[idx]
es = facesToEdges(tm.faces[fIdxs[fIdxs >= 0]])
p = edgesToPath(es[np.all(es!=idx,axis=1)])
es = np.transpose([p, np.roll(p, -1)])
es = np.vstack([es, es[:,::-1]])
es = np.vstack([es, np.transpose([p, [idx]*len(p)])])
self.nIdxs[idx,:,:len(es)] = es.T
em = bLoop[es].all(axis=1)
em[:len(p)*2] = True
cM = igl.cotmatrix(self.verts[col.indices], reIndexIndices(tm.faces[fIdxs[fIdxs>=0]]))
fes = reIndexIndices(es)
self.nWeights[idx,:len(es)] = cM[fes[:,0], fes[:,1]]
self.nWeights[idx,:len(es)] *= (1+em)
self.Ralpha = self.Ralpha.reshape(self.nVerts, -1, 1, 1)
self.R = np.repeat([np.eye(self.nDim, dtype = np.float32)], self.nVerts, axis=0)
self.zAll = np.zeros((self.nVerts, self.nDim))
self.uAll = np.zeros((self.nVerts, self.nDim))
self.rhoAll = np.full(self.nVerts, self.cfgParams.rhoInit)
self.zsAll = np.zeros((self.nVerts, self.nDim, self.nDim)) + self.pStress
self.usAll = np.zeros((self.nVerts, self.nDim, self.nDim))
self.rhosAll = np.full((self.nVerts, self.nDim), self.cfgParams.rhoInit)
self.lambdasAll = np.ones_like(self.vMass) * self.cfgParams.pLambda
self.energyVec = np.zeros(self.nVerts, dtype = np.float32)
self.truHoods = (self.verts[self.nIdxs[:,0]] - self.verts[self.nIdxs[:,1]]) * self.nWeights[:,:,None]
def optimize(self, printProgress = True):
self.cVertss = list(self.cVertss)
pStressE = self.pStressE[:,0] / self.pStressE[:,0].mean()
w = pStressE ** 0.5
w = np.ones_like(pStressE)
w /= w.sum()
stressWeights = icdf(self.vmStress)
ctsStr = lambda cts: (', cts: ' + ', '.join(['%f']*4))%tuple([cts.min(), cts.max(), np.dot(w, cts), (cts < np.deg2rad(0.1)).sum()/self.nVerts])
start = time()
for i in range(abs(self.cfgParams.maxIters)):
self.lambdasAll[:] = self.cfgParams.lambdaFun(i)
if self.cfgParams.useNormals:
self.lambdasAll[self.aIdxs] *= 5
if not i:
heatScals = computeRotationAngles(innerNxM(self.pStress, self.R)) if self.nDim == 2 else computeMinEulerAngles(np.eye(3), innerNxM(self.pStress, self.R))
self.ctss = [heatScals]
if self.cfgParams.stressWeighted:
self.lambdasAll *= (1 + stressWeights ** 0.5)
st = time()
k = self.localStep()
self.tLocal.append(time() - st)
st = time()
self.globalStep()
self.tGlobal.append(time() - st)
self.relativeDeltaV = np.max(np.abs(self.cVertss[-1] - self.cVertss[-2])) / np.max(np.abs(self.cVertss[-1] - self.verts))
heatScals = computeRotationAngles(innerNxM(self.pStress, self.R)) if self.nDim == 2 else computeMinEulerAngles(np.eye(3), innerNxM(self.pStress, self.R))
self.ctss.append(heatScals)
if printProgress:
if not i:
print("It:% 4s (%3d), E: %f, rDV: %f, tL: %f, tG: %f"%("X", k, self.energy, self.relativeDeltaV, self.tLocal[-1], self.tGlobal[-1]) + ctsStr(self.ctss[0]))
print("It:% 4d (%3d), E: %f, rDV: %f, tL: %f, tG: %f"%(i, k, self.energy, self.relativeDeltaV, self.tLocal[-1], self.tGlobal[-1]) + ctsStr(self.ctss[-1]))
if (self.cfgParams.maxIters > 0 and self.relativeDeltaV < self.cfgParams.stopDelta):
break
tL = sum(self.tLocal)
tG = sum(self.tGlobal)
tS = tL+tG
print('Time (s, l, g): %f, %f, %f'%(tS, tL, tG))
print('Time (s, l, g): %f, %f, %f'%(tS/tS, tL/tS, tG/tS))
self.cVertss = np.float64(self.cVertss)
def globalStep(self):
R = self.Ralpha * np.take(self.R, self.nIdxs[:,0], axis=0) + (1-self.Ralpha) * np.take(self.R, self.nIdxs[:,1], axis=0)
rhs = np.sum((R @ self.truHoods[...,None]).squeeze(), axis=1)
self.cVerts = self.solver(self.L @ rhs + self.cRHS) if self.cfgParams.useConstraints else self.solver(self.L.T @ rhs)
if not self.cfgParams.useConstraints:
cVerts = self.cVerts - self.cVerts.min(axis=0)
self.cMass = igl.massmatrix(self.cVerts, self.ts, 0).diagonal()
cVerts *= (self.volumeInit / self.cMass.sum()) ** (1/self.nDim)
self.cVerts = cVerts + (self.centerInit - cVerts.max(axis=0)/2)
self.cVertss.append(self.cVerts.copy())
def localStep(self):
self.energyVec *= 0
rotHoods = (self.cVerts[self.nIdxs[:,0]] - self.cVerts[self.nIdxs[:,1]])
Ss = np.transpose(self.truHoods, axes=[0,2,1]) @ rotHoods
msk = np.ones(len(Ss), np.bool_)
zs = self.zsAll.copy()
us = self.usAll.copy()
ns = self.pStress.copy()
rhos = self.rhosAll.copy()
for k in range(self.cfgParams.maxIterADMM):
S = Ss[msk] + outer2dWeighted(ns[msk], zs[msk] - us[msk], rhos[msk]) #/ self.nDim
R = np.transpose(orthogonalizeOrientations(S), axes=[0,2,1])
# z step
zsOld = zs[msk] + 0
Rn = innerNxM(ns[msk], R)
zs[msk] = lassoShrink(Rn + us[msk], (self.lambdasAll[msk] * self.vMass[msk]).reshape(-1,1,1) / rhos[msk].reshape(-1,self.nDim,1))
# u step
us[msk] += Rn - zs[msk]
rNorms = norm(zs[msk] - Rn)
sNorms = norm(zs[msk] - zsOld) * rhos[msk]
# rho setup
rMsk = np.transpose([msk.copy()]*self.nDim)
rMsk[msk] *= rNorms > self.cfgParams.mu * sNorms
sMsk = np.transpose([msk.copy()]*self.nDim)
sMsk[msk] *= np.bitwise_and(np.bitwise_not(rMsk[msk]), sNorms > self.cfgParams.mu * rNorms)
rhos[rMsk] *= self.cfgParams.tao
us[rMsk] /= self.cfgParams.tao
rhos[sMsk] /= self.cfgParams.tao
us[sMsk] *= self.cfgParams.tao
# stopping
epsPri = np.sqrt(2.0 * self.nDim) * self.cfgParams.abstol + self.cfgParams.reltol * np.max([norm(Rn).max(axis=1), norm(zs[msk]).max(axis=1)], axis=0)
epsDual = np.sqrt(self.nDim) * self.cfgParams.abstol + self.cfgParams.reltol * norm(rhos[msk].reshape(-1,self.nDim,1) * us[msk]).max(axis=1)
bMs = np.bitwise_and(rNorms < epsPri.reshape(-1,1), sNorms < epsDual.reshape(-1,1)).all(axis=1)
bMsk = msk.copy()
bMsk[msk] *= bMs
if np.any(bMsk):
self.zsAll[bMsk] = zs[bMsk].copy()
self.usAll[bMsk] = us[bMsk].copy()
self.rhosAll[bMsk] = rhos[bMsk].copy()
self.R[bMsk] = R[bMs].copy()
RdVminusDU = innerNxM(self.truHoods[bMsk], R[bMs]) - rotHoods[bMsk]
self.energyVec[bMsk] = self.cfgParams.pLambda * self.vMass[bMsk] * np.sum(np.abs(Rn[bMs]), axis=2).max(axis=1)
#self.energyVec[bMsk] = self.lambdasAll[bMsk] * self.vMass[bMsk] * np.sum(np.abs(Rn[bMs]), axis=2).max(axis=1)
self.energyVec[bMsk] += (inner1d(RdVminusDU * self.nWeights[bMsk][...,None], RdVminusDU)**2).sum(axis=1) / 2
msk *= np.bitwise_not(bMsk)
if not msk.sum():
break
self.energy = np.sum(self.energyVec)
self.Rs.append(self.R.copy())
return k
def fillGrid(self):
print('Generating inner grid.')
m = self.cfgParams.gridRes
n = m + 1
e = self.cVerts.max(axis=0) - self.cVerts.min(axis=0)
ns = np.int32((e/e.min()) * n + 0.5)
if self.nDim == 2:
nx,ny = ns
rx = np.linspace(0, 1, nx)
ry = np.linspace(0, 1, ny)
verts = np.transpose([np.tile(rx, ny), np.repeat(ry, nx)])
else:
nx,ny,nz = ns
rx = np.linspace(0, 1, nx)
ry = np.linspace(0, 1, ny)
rz = np.linspace(0, 1, nz)
verts = np.transpose([np.tile(rx, ny*nz), np.repeat(np.tile(ry, nz), nx), np.repeat(rz, nx*ny)])
self.cDims = (e / ns)
offset = self.cDims * self.cfgParams.finalLayer
verts *= self.cVerts.max(axis=0) - self.cVerts.min(axis=0) - offset
verts += self.cVerts.min(axis=0) + offset/2
tVerts = self.cVerts[self.ts]
innerOnly = (self.cVerts, self.tris[self.boundaryTrisMask]) if self.nDim == 3 else None
bWeights, tIdxs = computeBarycentricWeights(verts, tVerts, self.nDim, innerOnly, self.cfgParams.hullTolerance * self.cDims.mean())
if self.nDim == 2:
vIdxs = lambda i: [i, i+1, i+1+nx, i+nx]
ijs = np.vstack([[i,j] for j in range(ny-1) for i in range(nx-1)])
cells = np.vstack([vIdxs(i+nx*j) for i,j in ijs])
mx, my = nx-1, ny-1
nIdxs = []
for x in [-1,0,1]:
for y in [-1,0,1]:
nIdxs.append(x + mx*y)
self.cnIdxs = np.arange(len(cells)).reshape(-1,1) + np.tile(nIdxs, len(cells)).reshape(-1,len(nIdxs))
else:
vIdxs = lambda i: [i, i+1, i+1+nx, i+nx, i+nx*ny, i+1+nx*ny, i+1+nx+nx*ny, i+nx+nx*ny]
ijks = np.vstack([[i,j,k] for k in range(nz-1) for j in range(ny-1) for i in range(nx-1)])
cells = np.vstack([vIdxs(i+nx*j+nx*ny*k) for i,j,k in ijks])
mx, my, mz = nx-1, ny-1, nz-1
nIdxs = []
for x in [-1,0,1]:
for y in [-1,0,1]:
for z in [-1,0,1]:
nIdxs.append(x + mx*y + mx*my*z)
self.cnIdxs = np.arange(len(cells)).reshape(-1,1) + np.tile(nIdxs, len(cells)).reshape(-1,len(nIdxs))
cMsk = np.all(tIdxs[cells] > -1, axis=1)
self.cMsk = cMsk.copy()
print('nCells:', self.cMsk.sum(), len(cells))
# filter hull intersects
if self.nDim == 2:
if not self.cfgParams.hullTolerance:
hullEdges = self.cVerts[self.edges[self.boundaryEdgeMask]]
for cIdx, (cM, cell) in enumerate(zip(cMsk, cells)):
if cM:
es = faceToEdges(cell)
for edge in faceToEdges(cell):
hMsk = edgesIntersect2D(hullEdges[:,0], hullEdges[:,1], verts[edge[0]], verts[edge[1]])
if hMsk.any():
cMsk[cIdx] = False
break
else:
if embreeFound and not self.cfgParams.hullTolerance:
embreeDevice = rtc.EmbreeDevice()
scene = rtcs.EmbreeScene(embreeDevice)
mesh = TriangleMesh(scene, self.cVerts[self.tris[self.boundaryTrisMask]])
for cIdx, (cM, cell) in enumerate(zip(cMsk, cells)):
if cM:
es = hexaToEdges(cell)
eDirs, eLens = normVec(verts[es[:,1]] - verts[es[:,0]], True)
dsts = scene.run(np.float32(verts)[es[:,0]], np.float32(eDirs), query='DISTANCE')
if any(dsts < eLens):
cMsk[cIdx] = False
# filter kernels
nStart = cMsk.sum()
validNeighborMsk = np.bitwise_and(self.cnIdxs >= 0, self.cnIdxs < len(cells))
i = 0
for i in range(len(cells)):
cMsks = cMsk[self.cnIdxs * validNeighborMsk] * validNeighborMsk
orphanMsk = np.bitwise_and(np.sum(cMsks, axis=1) == 1, cMsks[:,(self.nDim**3)//2])
if np.any(orphanMsk):
cMsk[orphanMsk] = False
continue
cOrder = np.argsort(np.sum(cMsks,axis=1))
for cIdx in cOrder:
if cMsk[cIdx]:
msk = cMsk[self.cnIdxs[cIdx] * validNeighborMsk[cIdx]] * validNeighborMsk[cIdx]
dps = np.dot(mKernels[self.nDim], msk)
cMsk[cIdx] = np.all(dps != 2)
if not cMsk[cIdx]:
break
else:
break
print('filtered:', nStart - cMsk.sum())
print('nCells:', cMsk.sum(), len(cells))
cells = cells[cMsk]
usedVertIdxs = np.unique(cells.ravel())
self.bWeights = bWeights[usedVertIdxs]
self.gVertTIdxs = tIdxs[usedVertIdxs]
self.gVerts = verts[usedVertIdxs]
self.cells = reIndexIndices(cells) # flip lr?
if self.nDim == 2:
self.cells = self.cells[:,::-1]
diffs = self.verts - self.cVerts
ts = self.ts[self.gVertTIdxs]
ds = diffs[ts]
t = innerVxM(self.bWeights, np.transpose(ds,axes=[0,2,1]))
self.gVertz = self.gVerts + t
self.gInnerMask = np.ones(len(t), np.bool_)
def computeHulls(self, cubed = False, updateHullIdxs = True):
if self.nDim == 2:
triHullEdges = filterForSingleEdges(facesToEdges(self.ts, False))
tLoops = edgesToPaths(triHullEdges)
qHullEdges = filterForSingleEdges(facesToEdges(self.cells, False))
qLoops = edgesToPaths(qHullEdges)
tVertIdxs = np.unique(flatten(tLoops))
qVertIdxs = np.unique(flatten(qLoops))
self.qLoops = qLoops
if updateHullIdxs:
if hasattr(self, 'hullElementHashs'):
self.innerHullElementHashs = self.hullElementHashs.copy()
self.innerHullVertIdxs = self.hullVertIdxs.copy()
self.hullElementHashs = cantorPiV(qHullEdges)
self.hullVertIdxs = qVertIdxs
tVerts = self.cVerts[tVertIdxs] if cubed else self.verts[tVertIdxs]
qVerts = self.gVerts[qVertIdxs] if cubed else self.gVertz[qVertIdxs]
tEdges = reIndexIndices(np.vstack([pathToEdges(tLoop) for tLoop in tLoops]))
qEdges = reIndexIndices(np.vstack([pathToEdges(qLoop) for qLoop in qLoops]))
return tVerts, tEdges, qVerts, qEdges
else:
tVerts = self.cVerts if cubed else self.verts
oTris = self.tris[self.boundaryTrisMask]
tVerts = tVerts[np.unique(oTris.ravel())]
tris = reIndexIndices(oTris)
quads = hexasToFaces(self.cells)
unq, inv, cnt = np.unique(cantorPiKV(quads), return_inverse = True, return_counts = True)
boundaryQuadsMask = cnt[inv] == 1
if updateHullIdxs:
if hasattr(self, 'hullElementHashs'):
self.innerHullElementHashs = self.hullElementHashs.copy()
self.innerHullVertIdxs = self.hullVertIdxs.copy()
self.hullElementHashs = cantorPiKV(quads[boundaryQuadsMask])
self.hullVertIdxs = np.unique(quads[boundaryQuadsMask].ravel())
qVerts = self.gVerts if cubed else self.gVertz
self.oQuads = quads[boundaryQuadsMask]
qVertIdxs = np.unique(self.oQuads.ravel())
qVerts = qVerts[np.unique(self.oQuads.ravel())]
quads = reIndexIndices(self.oQuads)
return tVerts, tris, qVerts, quads
def padHull(self):
print('Padding hull.')
if self.nDim == 2:
tVerts, tEdges, qVerts, qEdges = self.computeHulls(self.cfgParams.useTransform)
self.pVerts = PadObject.getDeformedVerts(tVerts, tEdges, qVerts, qEdges, 0.1, True)
newCells = []
qEs = np.vstack([pathToEdges(qLoop) for qLoop in self.qLoops])
pEdges = qEdges + len(self.gVertz)
for qe, pe in zip(qEs[:,::-1], pEdges):
newCells.append(np.concatenate([qe, pe]))
else:
tVerts, tris, qVerts, quads = self.computeHulls(self.cfgParams.useTransform)
self.pVerts = PadObject.getDeformedVerts(tVerts, tris, qVerts, quads[:,::-1], 0.1, self.cfgParams.useProjection)
newHullVertIdxs = np.unique(quads.ravel()) + len(self.gVertz)
newCells = []
for oq, pq in zip(self.oQuads, quads):
newCells.append(np.concatenate([oq[::-1], pq[::-1] + len(self.gVertz)]))
self.cells = np.vstack([self.cells, newCells])
bWeights, tIdxs = computeBarycentricWeights(self.pVerts, self.cVerts[self.ts], self.nDim)
if self.cfgParams.hullTolerance:
outMask = np.any(self.bWeights < 0, axis=1)
if self.nDim == 2:
self.gVerts[outMask] = self.cVerts[self.ts[self.gVertTIdxs[outMask]]].mean(axis=1)
self.gVertz[outMask] = self.verts[self.ts[self.gVertTIdxs[outMask]]].mean(axis=1)
else:
qNormals = igl.per_vertex_normals(self.pVerts, quadsToTris(np.vstack([quads[:,::-1], np.roll(quads[:,::-1], 1, axis=1)])))
self.gVerts[self.hullVertIdxs] = self.pVerts - qNormals * self.cDims.mean()/2
self.bWeights[self.hullVertIdxs] = bWeights
self.gVertTIdxs[self.hullVertIdxs] = tIdxs
diffs = self.verts - self.cVerts
ts = self.ts[self.gVertTIdxs]
ds = diffs[ts]
t = innerVxM(self.bWeights, np.transpose(ds,axes=[0,2,1]))
self.gVertz = self.gVerts + t
self.bWeights = np.vstack([self.bWeights, bWeights])
self.gVertTIdxs = np.concatenate([self.gVertTIdxs, tIdxs])
if self.cfgParams.useTransform:
self.gVerts = np.vstack([self.gVerts, self.pVerts])
diffs = self.verts - self.cVerts
ds = diffs[self.ts[tIdxs]]
t = innerVxM(bWeights, np.transpose(ds, axes=[0,2,1]))
self.pVertz = self.pVerts + t
else: # experimental
self.gVerts = np.vstack([self.gVerts, qVerts])
self.pVertz = self.pVerts
self.gVertz = np.vstack([self.gVertz, self.pVertz])
vms = self.vmStress[self.ts[self.gVertTIdxs]]
self.vmStrezz = inner1d(self.bWeights, vms)
self.gInnerMask = np.concatenate([self.gInnerMask, np.zeros(len(self.pVerts), np.bool_)])
self.computeHulls()
if self.cfgParams.hullSmoothing:
self.smoothenHull()
def smoothenHull(self):
es = facesToEdges(self.cells) if self.nDim == 2 else hexasToEdges(self.cells)
newVertPozs = np.zeros((len(self.innerHullVertIdxs), self.nDim), np.float32)
newVertPoss = np.zeros((len(self.innerHullVertIdxs), self.nDim), np.float32)
cells = self.cells if self.nDim == 2 else hexOrderBT2Dg(self.cells)
for k in range(30):
SJs = computeJacobians(self.gVertz[cells], True)
sjMask = SJs < (0.01 if type(self.cfgParams.hullSmoothing) is bool else self.cfgParams.hullSmoothing)
print('MSJ:', SJs.min(), sjMask.sum(), len(sjMask))
if not np.any(sjMask):
break
vIdxs = cells[sjMask]
vIdxz = np.unique(np.concatenate([es[np.any(es == vIdx, axis=1)].ravel() for vIdx in vIdxs]))
newVertPozs = np.zeros((len(vIdxz), self.nDim), np.float32)
newVertPoss = np.zeros((len(vIdxz), self.nDim), np.float32)
for i, vIdx in enumerate(vIdxz):
nIdxs = np.unique(es[np.any(es == vIdx, axis=1)].ravel())
ws = [0.75, 0.25]
if k >= 5 and vIdx in vIdxs:
optiPt = computeOptiPoint(vIdx, cells, self.gVertz)
if vIdx in self.hullVertIdxs:
ws = [0.5, 0.5]
else:
optiPt = self.gVertz[nIdxs].mean(axis=0)
newVertPozs[i] = np.dot(ws, [optiPt, self.gVertz[vIdx]])
newVertPoss[i] = self.gVerts[nIdxs].mean(axis=0)
self.gVertz[vIdxz] = newVertPozs
self.gVerts[vIdxz] = newVertPoss
# reproject outer hull
tVerts, tElems, qVerts, qElems = self.computeHulls(False, updateHullIdxs = False)
pVerts = PadObject.getDeformedVerts(tVerts, tElems, qVerts, qElems, 0.1, self.nDim < 3)
hMsk = np.bool_([hvi in vIdxz for hvi in self.hullVertIdxs])
self.gVertz[self.hullVertIdxs[hMsk]] = pVerts[hMsk]
else:
print('Untangling failed, MSJ:', SJs.min())
def evalScaledJacobians(self):
if not hasattr(self, 'cells'):
print('object has no cells yet')
return
cells = self.cells if self.nDim == 2 else hexOrderBT2Dg(self.cells)
self.SJs = computeJacobians(self.gVertz[cells], True)
return self.SJs
def evalAlignment(self):
self.gStress = np.zeros((len(self.gVertz), self.nDim, self.nDim), np.float32)
self.gStressE = np.zeros((len(self.gVertz), self.nDim), np.float32)
sMats = (self.sMats[self.ts[self.gVertTIdxs]] * self.bWeights.reshape(-1,self.nDim+1,1,1)).sum(axis=1)
for gIdx in range(len(self.gVertz)):
self.gStress[gIdx], self.gStressE[gIdx] = computePrincipalStress(sMats[gIdx])
if not self.isRestored:
self.gStress = innerNxM(self.gStress, [self.MinitRot])
self.gHeat = np.zeros(len(self.gStress), np.float32)
self.gNeighbors = []
gEdges = facesToEdges(self.cells) if self.nDim == 2 else hexasToEdges(self.cells)
for vIdx in range(len(self.gVerts)):
nEdges = gEdges[np.any(gEdges == vIdx, axis=1)]
self.gNeighbors.append(nEdges[nEdges != vIdx])
nDirs = normVec(self.gVertz[self.gNeighbors[-1]] - self.gVertz[vIdx])
self.gHeat[vIdx] = np.arccos(np.clip(np.abs(np.dot(self.gStress[vIdx], nDirs.T)),0,1).max(axis=0)).mean()
gHeat = np.rad2deg(self.gHeat)
for heat in [gHeat, gHeat[self.gInnerMask]]:
vals = [heat.min(), heat.max(), heat.mean()] + [(heat < d).sum()/len(heat) for d in [1,5,10]]
print('Ea: ' + ('%0.6f, '*len(vals))%tuple(vals))
def show(self, cubed = True, withField = True):
mlab.figure()
es = facesToEdges(self.ts) if self.nDim == 2 else tetsToEdges(self.ts)
eTris = toEdgeTris(es)
verts = self.cVerts if cubed else self.verts
verts = pad2Dto3D(verts) if self.nDim == 2 else verts
col = (1,1,1) if cubed else (0,0,0)
# BBox
M = (mat2Dto3D(self.MinitRot) if self.nDim == 2 else self.MinitRot) if cubed else np.eye(3)
vMin = np.dot(verts, M.T).min(axis=0)
vMax = np.dot(verts, M.T).max(axis=0)
vs = (cubeVerts/2 + 0.5) * (vMax-vMin) + vMin
vs = np.dot(vs, M)
bx,by,bz = vs.T
mlab.triangular_mesh(bx,by,bz, toEdgeTris(facesToEdges(sixCubeFaces)), representation = 'mesh', color = (0.5,0.5,0.5), tube_radius = self.sf*0.5)
# tris / tets
x,y,z = verts.T
mPlot = mlab.triangular_mesh(x, y, z, eTris, representation = 'mesh', color = col, tube_radius = self.sf/5)
if withField:
x,y,z = np.repeat(verts.T, self.nDim, axis=1)
pStress = np.vstack(innerNxM(self.pStress, self.R)) if cubed else self.pStress
sDirs = np.vstack(pStress)
u,v,w = pad2Dto3D(sDirs).T if self.nDim == 2 else sDirs.T
scals = np.tile(np.arange(self.nDim)[::-1], len(verts))
self.sPlot = mlab.quiver3d(x, y, z, u, v, w, scalars = scals)
self.sPlot.glyph.color_mode = 'color_by_scalar'
self.sPlot.glyph.glyph_source.glyph_source.glyph_type = 'dash'
self.sPlot.glyph.glyph_source.glyph_position = 'center'
# quads / hexas
if hasattr(self, 'gVertz'):
verts = pad2Dto3D(self.gVertz) if self.nDim == 2 else self.gVertz
x,y,z = np.repeat(verts.T, self.nDim, axis=1)
sDirs = np.vstack(self.gStress)
if not cubed and withField:
u,v,w = pad2Dto3D(sDirs).T if self.nDim == 2 else sDirs.T
scals = np.tile(np.arange(self.nDim), len(verts))
self.sPlot = mlab.quiver3d(x, y, z, u, v, w, scalars = scals)
self.sPlot.glyph.color_mode = 'color_by_scalar'
self.sPlot.glyph.glyph_source.glyph_source.glyph_type = 'dash'
self.sPlot.glyph.glyph_source.glyph_position = 'center'
self.sPlot.module_manager.scalar_lut_manager.lut.table = rgb2rgba(np.eye(3)*255)
if hasattr(self, 'cells'):
mSJs = computeJacobians(self.gVertz[self.cells], True)
self.gHeat[:] = 1
for c, mSJ in zip(self.cells, mSJs):
self.gHeat[c] = np.minimum(self.gHeat[c], mSJ)
eTris = toEdgeTris(facesToEdges(self.cells) if self.nDim == 2 else hexasToEdges(self.cells))
gVerts = self.gVerts if cubed else self.gVertz
x,y,z = pad2Dto3D(gVerts).T if self.nDim == 2 else gVerts.T
mPlot = mlab.triangular_mesh(x, y, z, eTris, scalars = self.gHeat, representation = 'mesh', tube_radius = self.sf)
def restoreInitShape(self):
if self.isRestored:
print('Skip Restore')
return
print('Restore ...')
queue = [self.verts, self.cVerts]
if hasattr(self, 'pVerts'):
queue += [self.pVerts, self.pVertz, self.gVerts, self.gVertz]
elif hasattr(self, 'gVerts'):
queue += [self.gVerts, self.gVertz]
for verts in queue:
verts[:] = np.dot(verts, self.MinitRot)
verts *= (self.BBinit[1]-self.BBinit[0]).max()/2
verts += self.BBinit.mean(axis=0)
if len(self.cVertss) > 1:
for cVerts in self.cVertss:
cVerts[:] = np.dot(cVerts, self.MinitRot)
self.cVertss *= (self.BBinit[1]-self.BBinit[0]).max()/2
self.cVertss[:] += self.BBinit.mean(axis=0)
self.pStress = innerNxM(self.pStress, [self.MinitRot.T])
if hasattr(self, 'gOri'):
self.gStress = innerNxM(self.gStress, [self.MinitRot.T])
self.gOri = innerNxM(self.gOri, [self.MinitRot.T])
self.sf = norm(self.verts.max(axis=0)-self.verts.min(axis=0)) / 1000
self.isRestored = True
def saveState(self):
np.savez_compressed(tmpDir + '%s%s.npz'%(self.cfgParams.objName, 'Cubed'), cVerts = self.cVerts, R = self.R, isRestored = self.isRestored, protocol = 2)
np.savez_compressed(tmpDir + '%s%s.npz'%(self.cfgParams.objName, 'State'), cVertss = self.cVertss, Rs = self.Rs, ctss = self.ctss, protocol = 2)
def loadState(self, light = False):
cubedName = tmpDir + '%s%s.npz'%(self.cfgParams.objName, 'Cubed')
if os.path.exists(cubedName):
fileIn = np.load(cubedName, allow_pickle=True, encoding = 'latin1')
self.cVerts = fileIn['cVerts']
self.R = fileIn['R']
self.isRestored = fileIn['isRestored']
stateName = tmpDir + '%s%s.npz'%(self.cfgParams.objName, 'State')
if not light and os.path.exists(stateName):
fileIn = np.load(stateName, allow_pickle=True, encoding = 'latin1')
self.cVertss = fileIn['cVertss']
self.Rs = fileIn['Rs']
self.ctss = fileIn['ctss']
def exportHexasToObj(self, fileName = None, cubed = False):
fileName = 'data/%sHexas.obj'%self.cfgParams.objName if fileName is None else fileName
verts = []
faces = []
gVerts = self.gVerts if cubed else self.gVertz
for i, cell in enumerate(self.cells):
verts.append(gVerts[np.unique(cell.ravel())])
faces.append(reIndexIndices(hexaToFaces(cell)) + i*8)
writeObjFile(fileName, np.vstack(verts), np.vstack(faces))
def exportHexasToPly(self, fileName = None, cubed = False):
fileName = 'data/%s%sHexas.ply'%(self.cfgParams.objName, 'Cubed' if cubed else 'Res') if fileName is None else fileName
verts = []
cols = []
faces = []
gVerts = self.gVerts if cubed else self.gVertz
for cIdx, cell in enumerate(self.cells):
vIdxs = cell[np.argsort(cell)]
verts.append(gVerts[vIdxs])
cols.append(self.gHeat[vIdxs])
faces.append(reIndexIndices(hexaToFaces(cell)) + cIdx * 8)
writePlyFile(fileName, np.vstack(verts), np.vstack(faces)[:,::-1], verticesColors = np.transpose([np.concatenate(cols)]*3))
def exportTetsToPly(self, fileName = None, cubed = False):
fileName = 'data/%s%sTets.ply'%(self.cfgParams.objName, 'Cubed' if cubed else 'Init') if fileName is None else fileName
vs = self.cVerts if cubed and hasattr(self, 'cVerts') else self.verts
verts = []
cols = []
tris = []
tTris = tetraToFaces(np.arange(4))
cs = self.ctss[-1] if cubed else self.ctss[0]
for tIdx, tet in enumerate(self.ts):
verts.append(vs[tet])
cols.append(cs[tet])
tris.append(tTris + tIdx * 4)
writePlyFile(fileName, np.vstack(verts), np.vstack(tris), verticesColors = np.transpose([np.concatenate(cols)]*3))
def exportHullsToObj(self, fileName = None, cubed = False):
fileNameTris = 'data/%sHullTris.obj'%self.cfgParams.objName if fileName is None else fileName
fileNameQuads = 'data/%sHullQuads.obj'%self.cfgParams.objName if fileName is None else fileName
tVerts, tris, qVerts, quads = self.computeHulls(cubed)
writeObjFile(fileNameTris, tVerts, tris)
writeObjFile(fileNameQuads, qVerts, quads)
def exportCellsToMesh(self, fileName = None):
fileName = 'data/%s%s'%(self.cfgParams.objName, '.ply' if self.nDim == 2 else '.mesh') if fileName is None else fileName
if self.nDim == 2:
writePlyFile(fileName, pad2Dto3D(self.gVertz), self.cells)
else:
writeMeshFile(fileName, self.gVertz, self.cells)
if hasattr(self, 'vmStrezz'):
np.savetxt(fileName.split('.')[0]+'.vms', self.vmStrezz)
if sys.platform == 'win32': # some TVTK issue on linux?
class CubeVisual(HasTraits):
p = Range(0, 1000, 0)
scene = Instance(MlabSceneModel, ())
meshPlot = Instance(PipelineBase)
view = View(Item('scene', editor=SceneEditor(scene_class=MayaviScene), show_label=False),
Group('p'), resizable=True)
def showPlot(self, cObj):
self.cObj = cObj
self.eTris = toEdgeTris(facesToEdges(self.cObj.ts) if self.cObj.nDim == 2 else tetsToEdges(self.cObj.ts))
self.configure_traits()
@on_trait_change('p, scene.activated')
def update_plot(self):
p = min(len(self.cObj.cVertss)-1, self.p)
cVerts = pad2Dto3D(self.cObj.cVertss[p]) if self.cObj.nDim == 2 else self.cObj.cVertss[p]
rStress = innerNxM(innerNxM(innerNxM(self.cObj.pStress, [self.cObj.MinitRot]), self.cObj.Rs[p]), [self.cObj.MinitRot.T])
R = np.vstack(rStress)
Rdirs = pad2Dto3D(R) if self.cObj.nDim == 2 else R
Rdirs *= self.cObj.pStressE.reshape(-1,1) ** 0.5
heatScals = self.cObj.ctss[p] if hasattr(self.cObj, 'ctss') else computeMinEulerAngles(np.eye(3), rStress)
M = self.cObj.MinitRot if self.cObj.isRestored else np.eye(self.cObj.nDim)
if self.cObj.nDim == 2:
vMin = np.dot(cVerts[:,:-1], M.T).min(axis=0)
vMax = np.dot(cVerts[:,:-1], M.T).max(axis=0)
boxVerts = (cubeVerts[:,:-1]/2 + 0.5) * (vMax-vMin) + vMin
boxVerts = pad2Dto3D(np.dot(boxVerts, M))
else:
vMin = np.dot(cVerts, M.T).min(axis=0)
vMax = np.dot(cVerts, M.T).max(axis=0)
boxVerts = (cubeVerts/2 + 0.5) * (vMax-vMin) + vMin
boxVerts = np.dot(boxVerts, M)
bx,by,bz = boxVerts.T
if not hasattr(self, 'vPlot'):
x,y,z = pad2Dto3D(self.cObj.verts).T if self.cObj.nDim == 2 else self.cObj.verts.T
self.vPlot = mlab.triangular_mesh(x, y, z, self.eTris, representation = 'mesh', color = (1,1,1), tube_radius = self.cObj.sf/5)
if hasattr(self.cObj, 'aIdxs'):
s = np.ones(len(self.cObj.aIdxs))
self.sPlot = mlab.quiver3d(x[self.cObj.aIdxs], y[self.cObj.aIdxs], z[self.cObj.aIdxs], s, s, s, color=(1,1,1), scale_factor = self.cObj.sf*2, mode='sphere')
self.sPlot.glyph.color_mode = 'color_by_scalar'
self.sPlot.glyph.glyph_source.glyph_position = 'center'
x,y,z = cVerts.T
self.uPlot = mlab.triangular_mesh(x, y, z, self.eTris, representation = 'mesh', color = (0,0,0), tube_radius = self.cObj.sf/5)
self.bPlot = mlab.triangular_mesh(bx, by, bz, toEdgeTris(facesToEdges(sixCubeFaces)), representation = 'mesh', color = (0.75,0.75,0.75), tube_radius = self.cObj.sf*0.5)
self.hPlot = mlab.quiver3d(x, y, z, x*0+1, y*0+1, z*0+1, scalars = heatScals, mode='sphere', scale_factor = self.cObj.sf*5)
self.hPlot.glyph.color_mode = 'color_by_scalar'
self.hPlot.glyph.glyph_source.glyph_position = 'center'
x,y,z = np.repeat(cVerts, self.cObj.nDim, axis=0).T
u,v,w = Rdirs.T
scals = np.tile(np.arange(self.cObj.nDim)[::-1], self.cObj.nVerts)
self.rPlot = mlab.quiver3d(x, y, z, u, v, w, scalars = scals)
self.rPlot.glyph.color_mode = 'color_by_scalar'
self.rPlot.glyph.glyph_source.glyph_source.glyph_type = 'dash'
self.rPlot.glyph.glyph_source.glyph_position = 'center'
cntr = np.repeat([self.cObj.centerInit], 3, axis=0)
x,y,z = pad2Dto3D(cntr).T if self.cObj.nDim == 2 else cntr.T
u,v,w = mat2Dto3D(self.cObj.MinitRot).T if self.cObj.nDim == 2 else self.cObj.MinitRot.T
self.mPlot = mlab.quiver3d(x, y, z, u, v, w, scalars = [3,2,1])
self.mPlot.glyph.color_mode = 'color_by_scalar'
self.mPlot.glyph.glyph_source.glyph_source.glyph_type = 'dash'
self.mPlot.glyph.glyph_source.glyph_position = 'center'
else:
self.uPlot.mlab_source.points = cVerts
self.bPlot.mlab_source.points = boxVerts
self.rPlot.mlab_source.points = np.repeat(cVerts, self.cObj.nDim, axis=0)
self.rPlot.mlab_source.vectors = Rdirs
self.hPlot.mlab_source.points = cVerts
self.hPlot.mlab_source.scalars = heatScals
if hasattr(self.cObj, 'aIdxs'):
self.sPlot.mlab_source.points = cVerts[self.cObj.aIdxs]
if __name__ == "__main__":
# create own config
cfgParams = ConfigParams(objName = 'femur',
pLambda = 2,
lambdaRamp = 100,
autoScale = True,
optiAlign = True,
useNormals = True,
stressWeighted = False,
useConstraints = False,
maxIters = 1000,
stopDelta = 0.0001,
gridRes = 14,
finalLayer = False,
useProjection = False,
useTransform = True,
hullTolerance = False,
hullSmoothing = True)
# or load an existing one
#cfgParams = ConfigParams.fromFile('data/bar2D.cfg')
#cfgParams = ConfigParams.fromFile('data/bunny2D.cfg')
#cfgParams = ConfigParams.fromFile('data/spot2D.cfg')
#cfgParams = ConfigParams.fromFile('data/cube.cfg')
cfgParams = ConfigParams.fromFile('data/femur.cfg')
#cfgParams = ConfigParams.fromFile('data/fertility.cfg')
#cfgParams = ConfigParams.fromFile('data/kitten.cfg')
#cfgParams = ConfigParams.fromFile('data/spot.cfg')
#cfgParams = ConfigParams.fromFile('data/venus.cfg')
#cfgParams = ConfigParams.fromFile('data/JEB.cfg')
#cfgParams = ConfigParams.fromFile('data/buddhaDown.cfg')
#cfgParams = ConfigParams.fromFile('data/buddhaBack.cfg')
#cfgParams = ConfigParams.fromFile('data/buddhaBelly.cfg')
# load the stress file
cObj = CubeObject(cfgParams)
# either
if True: # compute field deformation and save it
cObj.optimize()
cObj.saveState()
else: # or load it from previous run
cObj.loadState()
# compute inner grid and hull
cObj.fillGrid()
#cObj.exportHullsToObj()
cObj.padHull()
cObj.restoreInitShape()
cObj.evalAlignment()
SJs = cObj.evalScaledJacobians()
print('#hv:\t%d, #he:\t%d'%(len(cObj.gVertz), len(cObj.cells)))
print('SJs:', SJs.min(), SJs.max(), SJs.mean())
# save results as mesh
cObj.exportCellsToMesh()
# show results in mayavi
#cObj.show(True, False)
cObj.show(False, False)
