from femUtil import * # FEM constants youngs_modulus = 1000.0 poisson_ratio = 0.3 pLambda, pMu = lame_parameters(youngs_modulus, poisson_ratio) C = linear_stress(pLambda, pMu) class FemObject: def __init__(self, inputFile = None, loadResult = False): if inputFile is None: return self.materialized = False self.inputFile = inputFile self.fixFuns = [] self.flxFuns = [] self.frcVecs = [] self.funData = [] if '.stress' in inputFile: self.fromStressFile(inputFile) return with open(self.inputFile, 'r') as fh: for line in fh.readlines(): line = line.strip().split(' ') if '#' in line: continue if line[0] == 'file': self.meshFile = os.path.dirname(self.inputFile) + '/' + line[1] continue if line[0] == 'type': self.mType = line[1] continue if line[1] == 'd': d = int(line[2]) if '-' in line[2]: fData = (-d, float(line[-1])) f = lambda vs, fIdx, tol = eps: vs[:,self.funData[fIdx][0]] <= (self.funData[fIdx][1] + tol) else: fData = (d, float(line[-1])) f = lambda vs, fIdx, tol = eps: vs[:,self.funData[fIdx][0]] >= (self.funData[fIdx][1] - tol) elif line[1] == 'r': fData = (np.float32(line[2:-1]), float(line[-1])) f = lambda vs, fIdx, tol = eps: norm(vs - self.funData[fIdx][0]) < self.funData[fIdx][1] + tol if line[0] == 'fix': self.fixFuns.append(f) self.funData.append(fData) elif line[0] == 'flx': self.flxFuns.append(f) self.funData.append(fData) elif line[0] == 'vec': self.frcVecs.append(np.float32(line[1:])) self.frcVecs = np.float32(self.frcVecs) self.nDim = 2 if self.mType in ['tri', 'quad'] else 3 if loadResult: if self.nDim == 2: self.meshFile = self.inputFile.replace('.frc', '.ply') self.mType = 'quad' if self.nDim == 3: self.meshFile = self.inputFile.replace('.frc', '.mesh') self.mType = 'hex' if not os.path.exists(self.meshFile): print('Result does not yet exist.') vmsPath = self.meshFile.split('.')[0] + '.vms' if os.path.exists(self.meshFile.split('.')[0] + '.vms'): self.vmStress = np.loadtxt(vmsPath) self.mTypeInit = self.mType if self.mType in ['tri', 'quad']: if '.obj' in self.meshFile: vs, self.ts = loadObjFile(self.meshFile) else: vs, self.ts = loadPlyFile(self.meshFile) self.verts = vs[:,:-1] if self.mType == 'tri': self.meshInit = fem.MeshTri(self.verts.T, self.ts.T) self.eType = fem.ElementTriP1() elif self.mType == 'quad': self.meshInit = fem.MeshQuad(self.verts.T, self.ts.T) self.eType = fem.ElementQuad1() elif self.mType == 'tet': self.meshInit = fem.MeshTet.load(self.meshFile) self.eType = fem.ElementTetP1() elif self.mType == 'hex': self.meshInit = fem.MeshHex.load(self.meshFile) self.eType = fem.ElementHex1() self.setupFEM() def fromStressFile(self, filePath): self.verts, self.ts, self.vIdxsFlx, self.forceVecs, self.vIdxsFix, self.vmStress, self.pStress, self.pStressE, self.sMats = loadStressFile(filePath) self.nDim = self.verts.shape[1] if self.nDim == 2: if self.ts.shape[1] == 3: self.mType = 'tri' self.meshInit = fem.MeshTri(self.verts.T, self.ts.T) self.eType = fem.ElementTriP1() else: self.mType = 'quad' self.meshInit = fem.MeshQuad(self.verts.T, self.ts.T) self.eType = fem.ElementQuad1() else: if self.ts.shape[1] == 4: self.mType = 'tet' self.meshInit = fem.MeshTet(self.verts.T, self.ts.T) self.eType = fem.ElementTetP1() else: self.mType = 'hex' self.meshInit = fem.MeshHex(self.verts.T, self.ts.T) self.eType = fem.ElementHex1() self.mTypeInit = self.mType self.setupFEM() fVecs = np.zeros((self.nVerts, self.nDim), np.float32) fVecs[self.vIdxsFlx] = self.forceVecs self.f = fVecs.ravel() def setupFEM(self): self.sf = norm(self.meshInit.doflocs.max(axis=1)-self.meshInit.doflocs.min(axis=1)) / 1000 self.nVerts = len(self.meshInit.doflocs.T) self.nElems = len(self.meshInit.t.T) print("#v % 6d, #e % 6d"%(self.nVerts, self.nElems)) vs = self.meshInit.doflocs.T es = self.meshInit.edges.T if self.nDim == 3 else facesToEdges(self.meshInit.t.T) self.mel = norm(vs[es[:,0]] - vs[es[:,1]]).mean() self.eVector = fem.ElementVector(self.eType) self.basis = fem.Basis(self.meshInit, self.eVector, fem.MappingIsoparametric(self.meshInit, self.eType), 0) self.K = fem.asm(linear_elasticity(pLambda, pMu), self.basis) if self.nDim == 2: tris = self.ts if self.mType == 'tri' else facesToTris(self.ts) self.volInit = computeTriangleAreas(self.verts[tris], False).sum() else: verts = self.meshInit.doflocs.T if self.mType == 'tet': tets = self.meshInit.t.T else: hexas = self.meshInit.t.T tets = np.vstack([hexaToTetras(hexa) for hexa in hexOrderDg2BT(hexas)]) self.volInit = computeTetraVolumes(verts[tets]).sum() self.vol = self.volInit self.setupBoundaryMasks() self.setupBoundaryWeights() def setupBoundaryWeights(self): self.bWeights = np.zeros(self.nVerts, np.float32) bElements = self.meshInit.facets.T[self.meshInit.boundary_facets()] bElements = bElements[np.any(self.boundaryVertMask[bElements], axis=1)] if self.nDim == 2: bElementWeights = norm(self.verts[bElements[:,0]] - self.verts[bElements[:,1]]) / 2 elif self.mType == 'tet': bElementWeights = computeTriangleAreas(self.meshInit.doflocs.T[bElements]) / 3 elif self.mType == 'hex': bElementWeights = computeTriangleAreas(self.meshInit.doflocs.T[tetrasToFaces(bElements)]).reshape(-1,4).sum(axis=1) / 8 for eIdx, bElement in tqdm(enumerate(bElements), total = len(bElements), ascii=True, desc='boundaryWeights'): self.bWeights[bElement] += bElementWeights[eIdx] def materialize(self, alpha = 0.5, walled = False): self.materialized = [alpha, walled] if self.mType in ['tri', 'quad']: self.materializeTriQuad(alpha) elif self.mType == 'tet': self.materializeTet(alpha, walled) elif self.mType == 'hex': self.materializeHex(alpha, walled) def materializeOpti(self, alpha = 0.5, walled = False, weights = None): self.materialized = [alpha, walled] alpha = abs(alpha) ws = (self.vmStress/self.vmStress.max()) ** 2 if weights is None else weights if self.mType in ['tri', 'quad']: vertsOrig = self.verts facesOrig = self.ts edgesOrig = facesToEdges(facesOrig) self.materializeTriQuad(alpha, False, False) taus = optimizeTaus(self.volInit * alpha, ws, vertsOrig, [edgesOrig, facesOrig], self.ts, self.mTypeInit) self.verts = materializeFaceVerts(taus, ws, vertsOrig, edgesOrig, facesOrig, self.mTypeInit) elif self.mType == 'tet': vertsOrig = self.meshInit.doflocs.T edgesOrig = self.meshInit.edges.T trisOrig = self.meshInit.facets.T tetsOrig = self.meshInit.t.T self.materializeTet(alpha, walled, False, False) taus = optimizeTaus(self.volInit * alpha, ws, vertsOrig, [edgesOrig, trisOrig, tetsOrig], self.ts, self.mTypeInit, walled) self.verts = materializeTetVerts(taus, ws, vertsOrig, edgesOrig, trisOrig, tetsOrig, self.mTypeInit, walled) elif self.mType == 'hex': vertsOrig = self.meshInit.doflocs.T edgesOrig = self.meshInit.edges.T quadsOrig = self.meshInit.facets.T hexasOrig = self.meshInit.t.T self.materializeHex(alpha, walled, False, False) taus = optimizeTaus(self.volInit * alpha, ws, vertsOrig, [edgesOrig, quadsOrig, hexasOrig], self.ts, self.mTypeInit, walled) self.verts = materializeHexVerts(taus, ws, vertsOrig, edgesOrig, quadsOrig, hexasOrig, self.mTypeInit, walled) self.initMaterializedFEM() def initMaterializedFEM(self): if self.mTypeInit in ['tri', 'quad', 'cell']: self.mType = 'quad' self.meshInit = fem.MeshQuad(self.verts.T, self.ts.T) self.eType = fem.ElementQuad1() self.vol = computeTriangleAreas(self.verts[facesToTris(self.ts)], False).sum() else: self.mType = 'hex' self.meshInit = fem.MeshHex(self.verts.T, self.ts.T) self.eType = fem.ElementHex1() hexaTets = np.vstack([hexaToTetras(hexa) for hexa in hexOrderDg2BT(self.ts)]) self.vol = computeTetraVolumes(self.verts[hexaTets]).sum() self.SJs = self.evalScaledJacobians() if self.SJs.max() < 0: self.ts = self.ts[:,::-1] self.SJs *= -1 self.nVerts = len(self.meshInit.doflocs.T) self.nElems = len(self.meshInit.t.T) print('nV: %d, nE: %d'%(self.nVerts, self.nElems)) st = time() self.eVector = fem.ElementVector(self.eType) self.basis = fem.Basis(self.meshInit, self.eVector, fem.MappingIsoparametric(self.meshInit, self.eType), 0) print('basis', time() - st) st = time() self.K = fem.asm(linear_elasticity(pLambda, pMu), self.basis) print('K', time() - st) st = time() self.setupBoundaryWeights() print('weights', time() - st) def materializeTriQuad(self, alpha = 0.5, withVerts = True, withFEM = True): tris = self.ts if self.mType == 'tri' else facesToTris(self.ts) self.volInit = computeTriangleAreas(self.verts[tris], False).sum() verts = self.meshInit.doflocs.T faces = self.cs if self.mTypeInit == 'cell' else self.ts edges = facesToEdges(faces) if withVerts: if self.mType == 'tri': tau = (1 - np.sqrt(1 - alpha)) / 3 if self.mType == 'quad': tau = (1 - np.sqrt(1 - alpha)) / 2 self.verts = materializeFaceVerts(None, tau, verts, edges, faces, self.mType) boundaryVertMask = flatten([[self.boundaryVertMask[edge].all()] * 2 for edge in edges]) es = np.repeat(edges, 2, axis=0) es[1::2] = edges[:,::-1] edgeHashs = cantorPiV(es, False) quads = [] idxOffset = 0 for fIdx, face in enumerate(faces): m = len(face) for vIdx in range(m): heCL = cantorPiO(face[vIdx], face[(vIdx+1)%m]) heCR = cantorPiO(face[vIdx], face[(vIdx-1)%m]) vCLidx = np.where(edgeHashs == heCL)[0][0] + len(verts) vCRidx = np.where(edgeHashs == heCR)[0][0] + len(verts) vCIidx = idxOffset + vIdx + len(verts) + len(edges)*2 quads.append([face[vIdx], vCLidx, vCIidx, vCRidx]) # corner heLC = cantorPiO(face[(vIdx+1)%m], face[vIdx]) vLCidx = np.where(edgeHashs == heLC)[0][0] + len(verts) vLIidx = idxOffset + (vIdx+1)%m + len(verts) + len(edges)*2 quads.append([vCLidx, vLCidx, vLIidx, vCIidx]) # edge idxOffset += m boundaryVertMask += [False] * m self.ts = np.vstack(quads) self.boundaryVertMask = np.concatenate([self.boundaryVertMask, boundaryVertMask]) if withFEM: self.initMaterializedFEM() def materializeTet(self, alpha = 0.5, walled = False, withVerts = True, withFEM = True): def tetToEdges(tet): return np.int64(np.transpose([tet[[0,0,0,1,2,3]], tet[[1,2,3,2,3,1]]])) def tetToFaces(tet): return np.int64(tet[[[0,2,3], [0,3,1], [0,1,2], [1,3,2]]]) verts = self.meshInit.doflocs.T edges = self.meshInit.edges.T faces = self.meshInit.facets.T tets = self.meshInit.t.T self.volInit = computeTetraVolumes(verts[tets]).sum() if walled: tau = (1 - (1 - alpha)**(1/3)) / 4 numEdgeVerts = 3 numFaceVerts = 7 numCellVerts = 14 else: taus = {0.25: 0.0943680425, 0.5: 0.1448414, 0.75: 0.195888335} if alpha in taus.keys(): # exact values tau = taus[alpha] else: # approximation with magic numbers p = np.sqrt(42)/5 + 2 # 3.29615807 q = -np.pi/2 - 1 - 2**(1/3)/2 + np.sqrt(6)/2 + 3*np.sqrt(2)/2 # 0.14528211 tau = np.arcsin(alpha*2-1)/(p*np.pi)+q numEdgeVerts = 2 numFaceVerts = 3 numCellVerts = 4 if withVerts: self.verts = materializeTetVerts(None, tau, verts, edges, faces, tets, self.mType, walled) boundaryVertMask = [] eHashs = cantorPiV(np.int64(edges)) ePos = {} for i, (e, eHash) in enumerate(zip(edges, eHashs)): for j, vIdx in enumerate(e): ePos[(eHash, vIdx)] = numEdgeVerts*i + j if walled: ePos[(eHash, None)] = numEdgeVerts*i + 2 boundaryVertMask += [self.boundaryVertMask[e].all()] * numEdgeVerts fHashs = cantorPiKV(np.int64(faces)) fPos = {} for i, (fc, fHash) in enumerate(zip(faces, fHashs)): for j, vIdx in enumerate(fc): fPos[(fHash, (vIdx+1))] = numFaceVerts*i + j if walled: fPos[(fHash, -(vIdx+1))] = numFaceVerts*i + j + 3 if walled: fPos[(fHash, None)] = numFaceVerts*i + 6 boundaryVertMask += [self.boundaryVertMask[fc].all()] * numFaceVerts newHexas = [] for tIdx, vIdxs in tqdm(enumerate(tets), total = tets.shape[0], ascii=True, desc='materializing'): tEdges = tetToEdges(vIdxs) tEdgesHashs = cantorPiV(tEdges) eIdxs = [] for edge, eHash in zip(tEdges, tEdgesHashs): eIdxs += [ePos[(eHash, vIdx)] for vIdx in edge] if walled: eIdxs += [ePos[(eHash, None)]] eIdxs = np.int32(eIdxs) + len(verts) tFaces = tetToFaces(vIdxs) tFacesHashs = cantorPiKV(tFaces) fIdxs = [] for face, fHash in zip(tFaces, tFacesHashs): fIdxs += [fPos[(fHash, (vIdx+1))] for vIdx in face] if walled: fIdxs += [fPos[(fHash, -(vIdx+1))] for vIdx in face] fIdxs += [fPos[(fHash, None)]] fIdxs = np.int32(fIdxs) + len(verts) + len(edges) * numEdgeVerts cIdxs = np.arange(numCellVerts) + tIdx*numCellVerts + len(verts) + len(edges)*numEdgeVerts + len(faces)*numFaceVerts # corners if walled: newHexas += [[vIdxs[0], eIdxs[0], eIdxs[6], eIdxs[3], fIdxs[7], fIdxs[14], fIdxs[0], cIdxs[0]]] newHexas += [[eIdxs[1], vIdxs[1], fIdxs[9], fIdxs[15], eIdxs[16], eIdxs[9], cIdxs[1], fIdxs[21]]] newHexas += [[eIdxs[4], fIdxs[16], fIdxs[1], vIdxs[2], cIdxs[2], eIdxs[10], eIdxs[12], fIdxs[23]]] newHexas += [[fIdxs[2], cIdxs[3], eIdxs[7], eIdxs[13], fIdxs[8], fIdxs[22], vIdxs[3], eIdxs[15]]] else: newHexas += [[vIdxs[0], eIdxs[0], eIdxs[4], eIdxs[2], fIdxs[3], fIdxs[6], fIdxs[0], cIdxs[0]]] newHexas += [[eIdxs[1], vIdxs[1], fIdxs[5], fIdxs[7], eIdxs[11], eIdxs[6], cIdxs[1], fIdxs[9]]] newHexas += [[eIdxs[3], fIdxs[8], fIdxs[1], vIdxs[2], cIdxs[2], eIdxs[7], eIdxs[8], fIdxs[11]]] newHexas += [[fIdxs[2], cIdxs[3], eIdxs[5], eIdxs[9], fIdxs[4], fIdxs[10], vIdxs[3], eIdxs[10]]] # edges if walled: newHexas += [[eIdxs[0], eIdxs[2], fIdxs[7], fIdxs[14], fIdxs[11], fIdxs[19], cIdxs[0], cIdxs[4]]] newHexas += [[eIdxs[2], eIdxs[1], fIdxs[11], fIdxs[19], fIdxs[9], fIdxs[15], cIdxs[4], cIdxs[1]]] newHexas += [[eIdxs[3], fIdxs[14], fIdxs[0], eIdxs[5], cIdxs[0], fIdxs[18], fIdxs[5], cIdxs[6]]] newHexas += [[eIdxs[5], fIdxs[18], fIdxs[5], eIdxs[4], cIdxs[6], fIdxs[16], fIdxs[1], cIdxs[2]]] newHexas += [[eIdxs[6], fIdxs[7], eIdxs[8], fIdxs[0], fIdxs[12], cIdxs[0], fIdxs[4], cIdxs[7]]] newHexas += [[eIdxs[8], fIdxs[12], eIdxs[7], fIdxs[4], fIdxs[8], cIdxs[7], fIdxs[2], cIdxs[3]]] newHexas += [[fIdxs[1], cIdxs[2], fIdxs[3], eIdxs[12], cIdxs[9], fIdxs[23], eIdxs[14], fIdxs[24]]] newHexas += [[fIdxs[3], cIdxs[9], fIdxs[2], eIdxs[14], cIdxs[3], fIdxs[24], eIdxs[13], fIdxs[22]]] newHexas += [[fIdxs[16], fIdxs[17], cIdxs[2], eIdxs[10], cIdxs[5], eIdxs[11], fIdxs[23], fIdxs[25]]] newHexas += [[fIdxs[17], fIdxs[15], cIdxs[5], eIdxs[11], cIdxs[1], eIdxs[9], fIdxs[25], fIdxs[21]]] newHexas += [[fIdxs[9], eIdxs[16], fIdxs[10], cIdxs[1], eIdxs[17], fIdxs[21], cIdxs[8], fIdxs[26]]] newHexas += [[fIdxs[10], eIdxs[17], fIdxs[8], cIdxs[8], eIdxs[15], fIdxs[26], cIdxs[3], fIdxs[22]]] else: newHexas += [[eIdxs[0], eIdxs[1], fIdxs[3], fIdxs[6], fIdxs[5], fIdxs[7], cIdxs[0], cIdxs[1]]] newHexas += [[eIdxs[2], fIdxs[6], fIdxs[0], eIdxs[3], cIdxs[0], fIdxs[8], fIdxs[1], cIdxs[2]]] newHexas += [[eIdxs[4], fIdxs[3], eIdxs[5], fIdxs[0], fIdxs[4], cIdxs[0], fIdxs[2], cIdxs[3]]] newHexas += [[fIdxs[1], cIdxs[2], fIdxs[2], eIdxs[8], cIdxs[3], fIdxs[11], eIdxs[9], fIdxs[10]]] newHexas += [[fIdxs[8], fIdxs[7], cIdxs[2], eIdxs[7], cIdxs[1], eIdxs[6], fIdxs[11], fIdxs[9]]] newHexas += [[fIdxs[5], eIdxs[11], fIdxs[4], cIdxs[1], eIdxs[10], fIdxs[9], cIdxs[3], fIdxs[10]]] # faces if walled: newHexas += [[fIdxs[6], cIdxs[11], fIdxs[5], fIdxs[4], cIdxs[6], cIdxs[7], fIdxs[0], cIdxs[0]]] newHexas += [[fIdxs[6], cIdxs[11], fIdxs[4], fIdxs[3], cIdxs[7], cIdxs[9], fIdxs[2], cIdxs[3]]] newHexas += [[fIdxs[6], cIdxs[11], fIdxs[3], fIdxs[5], cIdxs[9], cIdxs[6], fIdxs[1], cIdxs[2]]] newHexas += [[fIdxs[13], cIdxs[12], fIdxs[12], fIdxs[11], cIdxs[7], cIdxs[4], fIdxs[7], cIdxs[0]]] newHexas += [[fIdxs[13], cIdxs[12], fIdxs[11], fIdxs[10], cIdxs[4], cIdxs[8], fIdxs[9], cIdxs[1]]] newHexas += [[fIdxs[13], cIdxs[12], fIdxs[10], fIdxs[12], cIdxs[8], cIdxs[7], fIdxs[8], cIdxs[3]]] newHexas += [[fIdxs[20], cIdxs[13], fIdxs[19], fIdxs[18], cIdxs[4], cIdxs[6], fIdxs[14], cIdxs[0]]] newHexas += [[fIdxs[20], cIdxs[13], fIdxs[18], fIdxs[17], cIdxs[6], cIdxs[5], fIdxs[16], cIdxs[2]]] newHexas += [[fIdxs[20], cIdxs[13], fIdxs[17], fIdxs[19], cIdxs[5], cIdxs[4], fIdxs[15], cIdxs[1]]] newHexas += [[fIdxs[27], cIdxs[10], fIdxs[26], fIdxs[25], cIdxs[8], cIdxs[5], fIdxs[21], cIdxs[1]]] newHexas += [[fIdxs[27], cIdxs[10], fIdxs[25], fIdxs[24], cIdxs[5], cIdxs[9], fIdxs[23], cIdxs[2]]] newHexas += [[fIdxs[27], cIdxs[10], fIdxs[24], fIdxs[26], cIdxs[9], cIdxs[8], fIdxs[22], cIdxs[3]]] boundaryVertMask += [False] * numCellVerts self.ts = np.int32(newHexas)[:,::-1] self.boundaryVertMask = np.concatenate([self.boundaryVertMask, boundaryVertMask]) if withFEM: self.initMaterializedFEM() def materializeHex(self, alpha = 0.5, walled = False, withVerts = True, withFEM = True): def hexaToEdges(hexa): return np.int64(np.transpose([hexa[[0,0,0,1,2,1,3,2,3,4,5,6]], hexa[[1,2,3,4,4,5,5,6,6,7,7,7]]])) def hexaToFaces(hexa): return np.int64(hexa[[[0,2,6,3], [0,3,5,1], [0,1,4,2], [6,7,5,3], [2,4,7,6], [1,5,7,4]]]) verts = self.meshInit.doflocs.T edges = self.meshInit.edges.T faces = self.meshInit.facets.T hexas = self.meshInit.t.T hexaTets = np.vstack([hexaToTetras(hexa) for hexa in hexOrderDg2BT(hexas)]) self.volInit = computeTetraVolumes(verts[hexaTets]).sum() if withVerts: if walled: tau = (1 - (1 - alpha)**(1/3)) / 2 else: tau = np.arcsin(alpha*2-1)/(2*np.pi)+1/4 self.verts = materializeHexVerts(None, tau, verts, edges, faces, hexas, self.mType, walled) boundaryVertMask = [] eHashs = cantorPiV(np.int64(edges)) ePos = {} for i, (e, eHash) in enumerate(zip(edges, eHashs)): for j, vIdx in enumerate(e): ePos[(eHash, vIdx)] = 2*i + j boundaryVertMask += [self.boundaryVertMask[e].all()] * 2 fHashs = cantorPiKV(np.int64(faces)) fPos = {} for i, (fc, fHash) in enumerate(zip(faces, fHashs)): for j, vIdx in enumerate(fc): fPos[(fHash, vIdx)] = 4*i + j boundaryVertMask += [self.boundaryVertMask[fc].all()] * 4 newHexas = [] for hIdx, vIdxs in tqdm(enumerate(hexas), total = hexas.shape[0], ascii=True, desc='materializing'): hEdges = hexaToEdges(vIdxs) hEdgesHashs = cantorPiV(hEdges) eIdxs = [] for edge, eHash in zip(hEdges, hEdgesHashs): eIdxs += [ePos[(eHash, vIdx)] for vIdx in edge] eIdxs = np.int32(eIdxs) + len(verts) hFaces = hexaToFaces(vIdxs) hFacesHashs = cantorPiKV(hFaces) fIdxs = [] for face, fHash in zip(hFaces, hFacesHashs): fIdxs += [fPos[(fHash, vIdx)] for vIdx in face] fIdxs = np.int32(fIdxs) + len(verts) + len(edges)*2 cIdxs = np.arange(8) + hIdx*8 + len(verts) + len(edges)*2 + len(faces)*4 # corners newHexas += [[vIdxs[0], eIdxs[0], eIdxs[2], eIdxs[4], fIdxs[8], fIdxs[4], fIdxs[0], cIdxs[0]]] newHexas += [[eIdxs[1], vIdxs[1], fIdxs[9], fIdxs[7], eIdxs[6], eIdxs[10], cIdxs[1], fIdxs[20]]] newHexas += [[eIdxs[3], fIdxs[11], vIdxs[2], fIdxs[1], eIdxs[8], cIdxs[2], eIdxs[14], fIdxs[16]]] newHexas += [[eIdxs[5], fIdxs[5], fIdxs[3], vIdxs[3], cIdxs[3], eIdxs[12], eIdxs[16], fIdxs[15]]] newHexas += [[fIdxs[10], eIdxs[7], eIdxs[9], cIdxs[4], vIdxs[4], fIdxs[23], fIdxs[17], eIdxs[18]]] newHexas += [[fIdxs[6], eIdxs[11], cIdxs[5], eIdxs[13], fIdxs[21], vIdxs[5], fIdxs[14], eIdxs[20]]] newHexas += [[fIdxs[2], cIdxs[6], eIdxs[15], eIdxs[17], fIdxs[19], fIdxs[12], vIdxs[6], eIdxs[22]]] newHexas += [[cIdxs[7], fIdxs[22], fIdxs[18], fIdxs[13], eIdxs[19], eIdxs[21], eIdxs[23], vIdxs[7]]] # edges newHexas += [[eIdxs[4], fIdxs[4], fIdxs[0], eIdxs[5], cIdxs[0], fIdxs[5], fIdxs[3], cIdxs[3]]] newHexas += [[eIdxs[2], fIdxs[8], eIdxs[3], fIdxs[0], fIdxs[11], cIdxs[0], fIdxs[1], cIdxs[2]]] newHexas += [[fIdxs[1], cIdxs[2], eIdxs[14], fIdxs[2], fIdxs[16], cIdxs[6], eIdxs[15], fIdxs[19]]] newHexas += [[fIdxs[3], cIdxs[3], fIdxs[2], eIdxs[16], cIdxs[6], fIdxs[15], eIdxs[17], fIdxs[12]]] newHexas += [[eIdxs[0], eIdxs[1], fIdxs[8], fIdxs[4], fIdxs[9], fIdxs[7], cIdxs[0], cIdxs[1]]] newHexas += [[fIdxs[11], fIdxs[10], eIdxs[8], cIdxs[2], eIdxs[9], cIdxs[4], fIdxs[16], fIdxs[17]]] newHexas += [[fIdxs[5], fIdxs[6], cIdxs[3], eIdxs[12], cIdxs[5], eIdxs[13], fIdxs[15], fIdxs[14]]] newHexas += [[cIdxs[6], cIdxs[7], fIdxs[19], fIdxs[12], fIdxs[18], fIdxs[13], eIdxs[22], eIdxs[23]]] newHexas += [[fIdxs[7], eIdxs[10], cIdxs[1], fIdxs[6], fIdxs[20], eIdxs[11], cIdxs[5], fIdxs[21]]] newHexas += [[fIdxs[9], eIdxs[6], fIdxs[10], cIdxs[1], eIdxs[7], fIdxs[20], cIdxs[4], fIdxs[23]]] newHexas += [[cIdxs[4], fIdxs[23], fIdxs[17], cIdxs[7], eIdxs[18], fIdxs[22], fIdxs[18], eIdxs[19]]] newHexas += [[cIdxs[5], fIdxs[21], cIdxs[7], fIdxs[14], fIdxs[22], eIdxs[20], fIdxs[13], eIdxs[21]]] # faces if walled: newHexas += [[fIdxs[0], cIdxs[0], fIdxs[1], fIdxs[3], cIdxs[2], cIdxs[3], fIdxs[2], cIdxs[6]]] newHexas += [[fIdxs[8], fIdxs[9], fIdxs[11], cIdxs[0], fIdxs[10], cIdxs[1], cIdxs[2], cIdxs[4]]] newHexas += [[fIdxs[4], fIdxs[7], cIdxs[0], fIdxs[5], cIdxs[1], fIdxs[6], cIdxs[3], cIdxs[5]]] newHexas += [[cIdxs[2], cIdxs[4], fIdxs[16], cIdxs[6], fIdxs[17], cIdxs[7], fIdxs[19], fIdxs[18]]] newHexas += [[cIdxs[3], cIdxs[5], cIdxs[6], fIdxs[15], cIdxs[7], fIdxs[14], fIdxs[12], fIdxs[13]]] newHexas += [[cIdxs[1], fIdxs[20], cIdxs[4], cIdxs[5], fIdxs[23], fIdxs[21], cIdxs[7], fIdxs[22]]] # inner block #newHexas += [[cIdxs[0], cIdxs[1], cIdxs[2], cIdxs[3], cIdxs[4], cIdxs[5], cIdxs[6], cIdxs[7]]] boundaryVertMask += [False] * 8 self.ts = np.int32(newHexas) self.boundaryVertMask = np.concatenate([self.boundaryVertMask, boundaryVertMask]) if withFEM: self.initMaterializedFEM() def setupBoundaryMasks(self): self.boundaryVertMask = np.zeros(self.nVerts, dtype=np.bool_) if self.mType in ['tri', 'quad']: es = facesToEdges(self.ts if hasattr(self, 'ts') else self.meshInit.t.T, False) unq, inv, cnt = np.unique(cantorPiKV(es), return_inverse = True, return_counts = True) self.boundaryEdgeMask = cnt[inv] == 1 self.boundaryVertMask[np.unique(es[self.boundaryEdgeMask].ravel())] = True if self.mType in ['hex', 'tet']: ts = self.ts if hasattr(self, 'ts') else self.meshInit.t.T fcs = tetrasToFaces(ts) if self.mType == 'tet' else hexasToFaces(ts) unq, inv, cnt = np.unique(cantorPiKV(fcs), return_inverse = True, return_counts = True) self.boundaryFaceMask = cnt[inv] == 1 self.boundaryVertMask[self.meshInit.boundary_nodes()] = True def computeCompliance(self, complianceOnly = False): fVecs = np.zeros((self.nVerts, self.nDim), np.float32) if not hasattr(self, 'forceVecs'): vs = self.meshInit.doflocs.T funIdx = 0 fixMsk = self.boundaryVertMask * 0 for fixFun in self.fixFuns: msk = fixFun(vs, funIdx, self.mel/1000) fixMsk = np.bitwise_or(np.bitwise_and(msk, self.boundaryVertMask), fixMsk) funIdx += 1 self.vIdxsFix = np.where(fixMsk)[0].tolist() flxMsk = self.boundaryVertMask * 0 flxIdxs = [] for fIdx, flxFun in enumerate(self.flxFuns): msk = flxFun(vs, funIdx, self.mel/1000) bMsk = np.bitwise_and(msk, self.boundaryVertMask) flxMsk = np.bitwise_or(bMsk, flxMsk) fvIdxs = np.where(bMsk)[0] flxIdxs.append(fvIdxs) w = self.bWeights[fvIdxs] #* 0 + 1 fVecs[fvIdxs] = self.frcVecs[fIdx] * w.reshape(-1,1) / w.sum() funIdx += 1 numVertsPerForce = np.int32(list(map(len, flxIdxs))) self.vIdxsFlx = np.where(flxMsk)[0].tolist() else: fVecs[self.vIdxsFlx] = self.forceVecs self.f = fVecs.ravel() print('#fix: %d, #flx: %d'%(len(self.vIdxsFix), len(self.vIdxsFlx))) if False: msk = np.ones((self.nVerts, self.nDim), np.bool_) msk[self.vIdxsFix] = False self.I = np.where(msk.ravel())[0] self.cnds = fem.condense(self.K, b = self.f, I = self.I) else: # same but a bit faster self.D = (np.tile(np.int32(self.vIdxsFix)*self.nDim,[self.nDim,1]).T + np.arange(self.nDim)).ravel() self.cnds = fem.condense(self.K, b = self.f, D = self.D) print('solving...') if pypardisoFound: self.u = fem.solve(*self.cnds, solver = solver_pypardiso) else: self.u = fem.solve(*self.cnds) print('done') self.compliance = np.dot(self.u, self.K.dot(self.u))/2 self.dCols = norm(self.u.reshape(-1,self.nDim)) if complianceOnly: return self.compliance # deformed mesh and stresses self.meshDeformed = self.meshInit.translated(self.u.reshape(-1,self.nDim).T) self.dg1 = self.basis.with_element(self.eType) self.u1 = self.basis.interpolate(self.u) C_sym_grad = C(sym_grad(self.u1)) self.C = np.transpose(C_sym_grad, axes=[2,3,0,1]) if pypardisoFound: s00 = dg1Project(self.dg1, C_sym_grad[0,0]) s01 = dg1Project(self.dg1, C_sym_grad[0,1]) s11 = dg1Project(self.dg1, C_sym_grad[1,1]) if self.mType in ['tet', 'hex']: s02 = dg1Project(self.dg1, C_sym_grad[0,2]) s12 = dg1Project(self.dg1, C_sym_grad[1,2]) s22 = dg1Project(self.dg1, C_sym_grad[2,2]) else: s00 = self.dg1.project(C_sym_grad[0,0]) s01 = self.dg1.project(C_sym_grad[0,1]) s11 = self.dg1.project(C_sym_grad[1,1]) if self.mType in ['tet', 'hex']: s02 = self.dg1.project(C_sym_grad[0,2]) s12 = self.dg1.project(C_sym_grad[1,2]) s22 = self.dg1.project(C_sym_grad[2,2]) self.sMats = np.zeros((self.nVerts,self.nDim,self.nDim), np.float32) cts = np.zeros(len(self.sMats)) for tIdx, t in enumerate(self.meshInit.t.T): self.sMats[t] += self.C[tIdx] cts[t] += 1 self.sMats /= cts.reshape(-1,1,1) dets = simpleDets(self.sMats) self.sMats *= simpleSign(dets.reshape(-1,1,1)) self.pStress = np.zeros_like(self.sMats) self.pStressE = np.zeros((len(self.pStress), self.nDim), np.float32) for i, sMat in enumerate(self.sMats): self.pStress[i], self.pStressE[i] = computePrincipalStress(sMat) if self.mType in ['tet', 'hex']: self.vmStress = np.sqrt(((s00-s11)**2 + (s11-s22)**2 + (s22-s00)**2 + 6*(s01**2+s12**2+s02**2))/2) else: self.vmStress = np.sqrt(s00**2 + s11**2 - s00*s11 + 3*s01**2) return self.compliance def evalScaledJacobians(self, save = False): verts = self.verts if hasattr(self, 'verts') else self.meshInit.doflocs.T cells = self.ts if hasattr(self, 'ts') else self.meshInit.t.T self.SJs = computeJacobians(verts[cells], True) if save: mSJs = self.SJs.min(axis=1) mSJs.sort() np.savetxt(self.getFileName() + '.msjs', np.transpose([np.linspace(0,1,len(mSJs)), mSJs]), fmt='%.6e') return self.SJs def showPs(self): # somehow only works if called before any use of mayavi if not polyscopeFound or self.nDim < 3: print('Either polyscope is not installed or model is 2D.') return vsInit = self.verts if hasattr(self, 'verts') else self.meshInit.doflocs.T if self.mType == 'tet': tetsInit = self.ts if hasattr(self, 'ts') else self.meshInit.t.T psMeshInit = ps.register_volume_mesh("input_tetInit", vsInit, tets = tetsInit, enabled=True) else: hexasInit = self.meshInit.t.T[:,[0,3,6,2,1,5,7,4]] psMeshInit = ps.register_volume_mesh("input_hexInit", vsInit, hexes = hexasInit, enabled=True) if hasattr(self, 'vmStress'): psMeshInit.add_scalar_quantity("vMises", self.vmStress, enabled=True) psMeshInit.add_scalar_quantity("icdf", icdf(self.vmStress), enabled=False) if hasattr(self, 'vStress'): psMeshInit.add_scalar_quantity("vStress", self.vStress, enabled=False) if hasattr(self, 'SJs'): hs = self.ts if hasattr(self, 'ts') else self.meshInit.t.T mSJs = np.ones(len(vsInit)) np.minimum.at(mSJs, hs.ravel(), self.SJs.ravel()) psMeshInit.add_scalar_quantity("mSJs", mSJs, enabled=False, cmap = 'spectral') if hasattr(self, 'meshDeformed'): vsTrans = self.meshDeformed.doflocs.T if self.mType == 'tet': psTrans = ps.register_volume_mesh("result_tetDeformed", vsTrans, tets = tetsInit) else: psTrans = ps.register_volume_mesh("result_hexDeformed", vsTrans, hexes = hexasInit) psTrans.add_scalar_quantity("vMovement", norm(self.u.reshape(-1,3))) psTrans.add_scalar_quantity("vMises", self.vmStress, enabled=True) psTrans.add_scalar_quantity("icdf", icdf(self.vmStress), enabled=False) psTrans.add_scalar_quantity("pSe0", normZeroToOne(self.pStressE[:,0]), enabled=False) psTrans.add_scalar_quantity("pSe1", normZeroToOne(self.pStressE[:,1]), enabled=False) psTrans.add_scalar_quantity("pSe2", normZeroToOne(self.pStressE[:,2]), enabled=False) ps.show() def show(self, showField = True, showVerts = False, showVectors = False, showDeformed = True, showLabels = False): if self.mType in ['tri', 'quad']: elementsInit = self.ts if hasattr(self, 'ts') else self.meshInit.t.T es = facesToEdges(elementsInit) elif self.mType == 'tet': elementsInit = self.ts if hasattr(self, 'ts') else self.meshInit.t.T es = tetsToEdges(elementsInit) elif self.mType == 'hex': elementsInit = self.meshInit.t.T[:,[0,3,6,2,1,5,7,4]] es = hexasToEdges(elementsInit) eTris = toEdgeTris(es) verts = self.verts if hasattr(self, 'verts') else self.meshInit.doflocs.T vsInit = pad2Dto3D(verts) if self.nDim == 2 else verts x,y,z = vsInit.T scals = icdf(self.vmStress) ** 0.5 if hasattr(self, 'vmStress') else np.zeros_like(x) # init wireframe if self.nDim == 2: self.mPlot = mlab.triangular_mesh(x, y, z, facesToTris(elementsInit), representation = 'surface', scalars = scals) self.mPlot = mlab.triangular_mesh(x, y, z, eTris, representation = 'mesh', tube_radius=self.sf*0.1, color = (0,0,0)) else: self.mPlot = mlab.triangular_mesh(x, y, z, eTris, representation = 'mesh', tube_radius=self.sf*2, scalars = scals) # deformed wireframe if showDeformed and hasattr(self, 'meshDeformed'): vsTrans = self.meshDeformed.doflocs.T if self.mType in ['tet', 'hex'] else pad2Dto3D(self.meshDeformed.doflocs.T) xt,yt,zt = vsTrans.T if self.nDim == 2: self.dPlot = mlab.triangular_mesh(xt, yt, zt, facesToTris(elementsInit), representation = 'surface', scalars = scals) else: self.dPlot = mlab.triangular_mesh(xt, yt, zt, eTris, representation = 'mesh', tube_radius = self.sf*2, scalars = scals) # verts if showVerts: s = np.ones_like(x) self.vPlot = mlab.quiver3d(x, y, z, s, s, s, scalars = scals, scale_factor = 0.05, mode='sphere') self.vPlot.glyph.color_mode = 'color_by_scalar' self.vPlot.glyph.glyph_source.glyph_position = 'center' # labels if showLabels: # very slow, impractical for large models mlab.gcf().scene.disable_render = True # should speedup rendering but has no effect somehow for vIdx, v in enumerate(vsInit): mlab.text3d(v[0], v[1], v[2], str(vIdx), scale = (self.sf*10,self.sf*10,self.sf*10)) mlab.gcf().scene.disable_render = False # deformation vectors if showVectors and hasattr(self, 'u'): uVecs = pad2Dto3D(self.u.reshape(-1,self.nDim)) u,v,w = uVecs.T s = norm(uVecs) self.uPlot = mlab.quiver3d(x, y, z, u, v, w, scalars = s, scale_factor = 1, mode='arrow') # stress field if showField and hasattr(self, 'pStress'): x,y,z = np.repeat(vsInit.T, self.nDim, axis=1) pStressScale = 1#np.clip(self.pStressE / np.percentile(self.pStressE, 99, axis=0), 0, 1).reshape(-1,3,1) ** 0.5 pStress = np.vstack(self.pStress * pStressScale) u,v,w = pStress.T if self.mType in ['tet', 'hex'] else pad2Dto3D(pStress).T scals = np.tile(np.arange(self.nDim)[::-1], len(vsInit)) 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' def showForces(self, showLabels = False): if self.mType in ['tri', 'quad']: elementsInit = self.meshInit.t.T es = facesToEdges(elementsInit) elif self.mType == 'tet': elementsInit = self.meshInit.t.T es = tetsToEdges(elementsInit) elif self.mType == 'hex': elementsInit = self.meshInit.t.T[:,[0,3,6,2,1,5,7,4]] es = hexasToEdges(elementsInit) eTris = toEdgeTris(es) vsInit = pad2Dto3D(self.meshInit.doflocs.T) if self.nDim == 2 else self.meshInit.doflocs.T x,y,z = vsInit.T # edges self.ePlot = mlab.triangular_mesh(x, y, z, eTris, representation = 'mesh', tube_radius=self.sf/5, color=(0,0,0)) # labels if showLabels: # very slow, impractical for large models mlab.gcf().scene.disable_render = True # should speedup rendering but has no effect somehow for vIdx, v in enumerate(vsInit): if self.boundaryVertMask[vIdx]: mlab.text3d(v[0], v[1], v[2], str(vIdx), scale = (self.sf*2.5,self.sf*2.5,self.sf*2.5)) mlab.gcf().scene.disable_render = False # vertices scals = np.zeros_like(x) scals[self.vIdxsFix] = 1 scals[self.vIdxsFlx] = 2 print('nFix: %d, nFlex: %d'%(len(self.vIdxsFix), len(self.vIdxsFlx))) s = np.ones_like(x) self.vPlot = mlab.quiver3d(x, y, z, s, s, s, scalars = scals, scale_factor = self.sf*2, mode='sphere') self.vPlot.glyph.color_mode = 'color_by_scalar' self.vPlot.glyph.glyph_source.glyph_position = 'center' # applied forces if hasattr(self, 'f'): x,y,z = vsInit[self.vIdxsFlx].T fVecs = pad2Dto3D(self.f.reshape(-1,self.nDim)) if self.nDim == 2 else self.f.reshape(-1,self.nDim) u,v,w = fVecs[self.vIdxsFlx].T self.fPlot = mlab.quiver3d(x, y, z, u, v, w, scale_factor = self.sf*50, mode='arrow') self.fPlot.glyph.glyph_source.glyph_position = 'tail' def exportToStress(self, fileName = None): fileName = self.getFileName(fileName).split('_')[0] + '.stress' verts = self.verts if hasattr(self, 'verts') else self.meshInit.doflocs.T cells = self.ts if hasattr(self, 'ts') else self.meshInit.t.T fixedIdxs = self.vIdxsFix forceIdxs = self.vIdxsFlx forceVecs = self.f.reshape(-1, self.nDim)[self.vIdxsFlx] sMats = np.transpose(self.sMats, axes=[1,2,0]) stress = sMats[[0,1,0],[0,1,1],:].T if self.nDim == 2 else sMats[[0,1,2,1,0,0],[0,1,2,2,2,1],:].T writeStressFile(fileName, verts, cells, forceIdxs, forceVecs, fixedIdxs, stress) def exportToMesh(self, fileName = None): fileName = self.getFileName(fileName) + ('.mesh' if self.nDim == 3 else '.obj') verts = self.verts if hasattr(self, 'verts') else self.meshInit.doflocs.T cells = self.ts if hasattr(self, 'ts') else self.meshInit.t.T if self.nDim == 2: writeObjFile(fileName, pad2Dto3D(verts), cells) elif self.mType == 'hex': writeMeshFile(fileName, verts, hexOrderDg2BT(cells[:,::-1])) def getFileName(self, fileName = None): if fileName is None: fileName = 'data/' + os.path.basename(self.inputFile).split('.')[0] fileName += '_' + self.mTypeInit if self.materialized: fileName += '_' + ('wall' if self.materialized[1] else 'beam') fileName += ('_%d'%(self.materialized[0]*100)) if self.materialized[0] > 0 else '_tau' return fileName if __name__ == "__main__": #inFile = 'data/bar2D.frc' #inFile = 'data/spot2D.frc' #inFile = 'data/bunny2D.frc' #inFile = 'data/cube.frc' inFile = 'data/femur.frc' #inFile = 'data/fertility.frc' #inFile = 'data/kitten.frc' #inFile = 'data/spot.frc' #inFile = 'data/venus.frc' #inFile = 'data/JEB.frc' #inFile = 'data/buddhaDown.frc' #inFile = 'data/buddhaBack.frc' #inFile = 'data/buddhaBelly.frc' fObj = FemObject(inFile, False) #fObj = FemObject(inFile, True) c0 = fObj.computeCompliance() print('c0:', c0) fObj.exportToStress() """ fObj.materialize(0.5, False) #fObj.materializeOpti(0.5, False, weights = icdf(fObj.vmStress)**2) c1 = fObj.computeCompliance() print('c1:', c1) print('c:', c0, c1, c1/c0) print('v:', fObj.volInit, fObj.vol, fObj.vol/fObj.volInit) """ #SJs = fObj.evalScaledJacobians() #print('msj:', fObj.SJs.min(), fObj.SJs.max(), SJs.mean()) fObj.showPs() fObj.show(False, showDeformed = False) #fObj.showForces()