# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE.txt', which is part of this source code package.
#
""" Class to treat Constants Rectilinear Finite Element Problems
"""
from typing import Optional
import numpy as np
from scipy.sparse import coo_matrix
import scipy.sparse.linalg as linalg
import scipy.linalg as denselinalg
import scipy.sparse as sps
from Muscat.Types import MuscatIndex, MuscatFloat, ArrayLike
import Muscat.Containers.ElementsDescription as ED
import Muscat.Containers.Mesh as Mesh
import Muscat.FE.FeaBase as FeaBase
from Muscat.LinAlg.MatVecOperations import DeleteRowsCols
from Muscat.Helpers.Logger import Info, Debug, Error
from Muscat.Containers.Filters.FilterObjects import NodeFilter
from Muscat.Containers.ConstantRectilinearMeshTools import GetMonoIndex
[docs]class ElementaryMatrix():
def __init__(self,dim=3, physics="disp"):
self.dim = dim
self.geoFactor = [1,1,1]
self.physics = physics
self.young = 1
self.poisson = 0.3
self.density = 1
self.thermalK = 1
[docs] def GetMassMatrix(self):
if self.physics == "disp":
from Muscat.FE.SymPhysics import MecaPhysics
phys = MecaPhysics(self.dim)
wform = phys.GetBulkMassFormulation()
return self.Integrate(["u_"+str(x) for x in range(self.dim)],wform)
elif self.physics == "thermal":
from Muscat.FE.SymPhysics import ThermalPhysics
wform = ThermalPhysics().GetBulkMassFormulation()
return self.Integrate(["t"],wform)
[docs] def GetTangentMatrix(self):
from Muscat.FE.SymPhysics import MecaPhysics,ThermalPhysics
if self.physics == "disp":
physics = MecaPhysics(dim=self.dim)
wform = physics.GetBulkFormulation(self.young,self.poisson)
elif self.physics == "thermal":
physics = ThermalPhysics()
physics.spaceDimension = self.dim
wform = physics.GetBulkFormulation(alpha=self.thermalK)
else:
raise(Exception("Physics not supported"))
return self.Integrate(physics.GetPrimalNames() ,wform)
[docs] def Integrate(self,primalNames,wform):
from Muscat.FE.FETools import GetElementaryMatrixForFormulation
if self.dim == 3:
el = ED.Hexahedron_8
else:
el = ED.Quadrangle_4
KGeneric,rhs = GetElementaryMatrixForFormulation(el,wform, unknownNames = primalNames, geoFactor= self.geoFactor)
n = len(primalNames)
permutation = np.array([x*ED.numberOfNodes[el]+y for y in range(ED.numberOfNodes[el]) for x in range(n)])
rhs = rhs[permutation]
KGeneric = KGeneric[permutation,:]
KGeneric = KGeneric[:,permutation]
return KGeneric
[docs]class Fea(FeaBase.FeaBase):
def __init__(self,spaceDim=3):
super().__init__(spaceDim=spaceDim)
self.linearSolver = "EigenCG"
self.writer = None
self.minthreshold = 0.9e-3
self.tol = 1.e-6
self.dofpernode = 1
self.init = False
self.dirichletBcs = []
self.neumannBcs = []
self.neumannNodalBcs = []
self.young = 1.
self.poisson = 0.3
self.density = 1.
[docs] def AddDirichletCondition(self, nodeFilter, mask, forceFunction):
self.dirichletBcs.append( (nodeFilter, mask, forceFunction) )
[docs] def AddNeumannCondition(self, nodeFilter, mask, forceFunction):
self.neumannBcs.append( (nodeFilter, mask, forceFunction) )
[docs] def AddNeumannNodalCondition(self, nodeFilter, mask, forceFunction):
self.neumannNodalBcs.append( (nodeFilter, mask, forceFunction) )
[docs] def BuildProblem(self, dofpernode = None, KOperator= None, MOperator= None, neumannNodalBcs= None):
if self.init == True:
return
self.init = True
if self.mesh.props.get("IsConstantRectilinear",False) == False :
raise Exception("Must be a ConstantRectilinear mesh type ") #pragma: no cover
self.outer_v = []
if dofpernode is not None:
self.dofpernode = dofpernode
self.nodesPerElement = 2**self.mesh.GetElementsDimensionality()
if KOperator is not None:
self.KE = KOperator
self.ME = MOperator
else:
self.myElem = ElementaryMatrix(dim = self.mesh.GetElementsDimensionality())
self.myElem.young = self.young
self.myElem.poisson = self.poisson
self.myElem.geoFactor = self.mesh.props.get("spacing")
self.KE = self.myElem.GetTangentMatrix()
self.ME = self.myElem.GetMassMatrix()
# dofs:
self.ndof = self.dofpernode * self.mesh.GetNumberOfNodes()
# FE: Build the index vectors for the for coo matrix format
Debug("Building Connectivity matrix")
nbBulkElements = self.mesh.GetNumberOfElements(self.mesh.GetElementsDimensionality())
self.edofMat = np.zeros((nbBulkElements, self.nodesPerElement*self.dofpernode), dtype=MuscatIndex)
Debug("Building Connectivity matrix 2")
bulkElements = self.mesh.GetElementsOfType([ED.Quadrangle_4,ED.Hexahedron_8][self.mesh.GetElementsDimensionality()-2])
bulkConnectivity = bulkElements.connectivity
for i in range(nbBulkElements):
conn = bulkConnectivity[i,:]
self.edofMat[i, :] = np.array([(conn*self.dofpernode+y) for y in range(self.dofpernode)]).ravel('F')
Debug("Building Connectivity matrix Done")
self.iK = None
self.jK = None
self.fixedValues = np.zeros((self.ndof, 1), dtype=MuscatFloat)
Debug("Treating Dirichlet 1/4")
self.fixed = np.zeros(self.ndof, dtype=bool)
for ef, mask, forceFunction in self.dirichletBcs:
indexs = ef.GetNodesIndices(mesh = self.mesh)
values = forceFunction(self.mesh.nodes[indexs,:])
for i, m in enumerate(mask):
if m:
self.fixed[indexs*self.dofpernode+i] = True
self.fixedValues[indexs*self.dofpernode+i] = values[i]
self.free = np.ones(self.ndof, dtype=bool)
self.free[self.fixed] = False
# Solution and RHS vectors
self.f = np.zeros((self.ndof, 1), dtype=MuscatFloat)
self.u = np.zeros((self.ndof, 1), dtype=MuscatFloat)
Debug("Treating Dirichlet Done")
Debug("Treating Neumann")
for ef, mask, forceFunction in self.neumannBcs:
nIndices = ef.GetNodesIndices(mesh = self.mesh)
if len(nIndices) == 0:
continue
z = np.zeros((self.mesh.GetNumberOfNodes(),))
f = np.zeros((self.ndof, 1), dtype=MuscatFloat)
values = forceFunction(self.mesh.nodes[nIndices,:])
for i, m in enumerate(mask):
if m:
z[nIndices] = 1.
f[nIndices*self.dofpernode + i] += values[i]
eff = np.clip((np.sum(z[bulkConnectivity],axis=1) ),0, 1)
MassMatrix = self.BuildMassMatrix(eff)
self.f[:] += MassMatrix*f[:]
if len (self.neumannNodalBcs) > 0:
nodal_f = np.zeros((self.ndof, 1), dtype=MuscatFloat)
for ef, mask, forceFunction in self.neumannNodalBcs:
nIndices = ef.GetNodesIndices(mesh = self.mesh)
values = forceFunction(self.mesh.nodes[nIndices,:])
for i, m in enumerate(mask):
if m:
nodal_f[nIndices*self.dofpernode + i] += values[i]
self.f[:,0] += nodal_f[:,0]
Debug("Treating Neumann Done")
self.eed = np.zeros(nbBulkElements)
[docs] def AssemblyMatrix(self,Op, Eeff = None):
nbBulkElements = self.mesh.GetNumberOfElements(self.mesh.GetElementsDimensionality())
if Eeff is None:
Eeff = np.ones(nbBulkElements)
sM = ((Op.flatten()[np.newaxis]).T * Eeff.ravel()).flatten(order='F')
self.GenerateIJs()
M = coo_matrix((sM, (self.iK, self.jK)), shape=(self.ndof, self.ndof),dtype=MuscatFloat)
else:
Debug(" Eeff is known")
bool_Eeff = (Eeff>=self.minthreshold)
nEeff = Eeff[bool_Eeff]
sM = ((Op.flatten()[np.newaxis]).T * nEeff.ravel()).flatten(order='F')
one = np.ones((self.nodesPerElement*self.dofpernode, 1), dtype=MuscatIndex)
local_iK = np.kron(self.edofMat[bool_Eeff,:], one).flatten()
one.shape = (1,self.nodesPerElement*self.dofpernode)
local_jK = np.kron(self.edofMat[bool_Eeff,:], one).flatten()
M = coo_matrix((sM, (local_iK, local_jK)), shape=(self.ndof, self.ndof),dtype=MuscatFloat)
return M.tocsr()
[docs] def BuildTangentMatrix(self, Eeff = None):
Debug("BuildTangentMatrix")
res = self.AssemblyMatrix(self.KE, Eeff)
Debug("BuildTangentMatrix Done")
return res
[docs] def BuildMassMatrix(self, Eeff = None):
Debug("BuildMassMatrix")
res = self.AssemblyMatrix(self.ME, Eeff)
Debug("BuildMassMatrix Done")
return res
[docs] def Solve(self, Eeff=None):
# hack to integrate complex boundary condition in the mesh
if hasattr(self,"mecaPhysics") and (self.mecaPhysics is not None):
from Muscat.FE.UnstructuredFeaSym import UnstructuredFeaSym
prob = UnstructuredFeaSym(spaceDim =self.spaceDim)
prob.physics.append(self.mecaPhysics)
prob.SetMesh(self.mesh)
self.mecaPhysics.ComputeDofNumberingFromConnectivity(self.mesh)
prob.ComputeDofNumbering(self.mesh)
k,f = prob.GetLinearProblem(computeK=False)
if f is not None:
f.shape = (self.spaceDim,len(f)//self.spaceDim)
f = f.ravel(order='F')
self.f[:,0] += f
for nf,data in self.mecaPhysics.extraRHSTerms:
nids = nf.GetNodesIndices(self.mesh)
size = 0
offset = []
for n in self.mecaPhysics.numberings:
offset.append(size)
size += n.size
for i,val in enumerate(data):
inumbering = self.mecaPhysics.numberings[i]
if val == 0 :
continue
for nid in nids:
dof = inumbering.GetDofOfPoint(nid)+offset[i]
self.f[dof,0] = val
self.mecaPhysics = None
Debug("Construction of the tangent matrix")
K = self.BuildTangentMatrix(Eeff)
zerosdof = np.where(K.diagonal()== 0 )[0]
Info("Number of active nodes : " + str(self.ndof-len(zerosdof) ) + " of " + str(self.ndof) + " "+ str(float(len(zerosdof)*100.)/self.ndof)+ "% of empty dofs" )
Kones = coo_matrix( (np.ones((len(zerosdof),) ) ,(zerosdof,zerosdof)), shape =(self.ndof, self.ndof)).tocsr()#(self.dofpernode,self.dofpernode))
K = (K.tocsr() + Kones.tocsr()).tocsr()
# Remove constrained dofs from matrix
Debug(" Delete fixed Dofs")
[K, rhsfixed] = DeleteRowsCols(K, self.fixed, self.fixed, self.fixedValues)
Debug(" Start solver (" + str(self.linearSolver) + ")")
rhs = self.f[self.free, 0]-rhsfixed[self.free, 0]
self.u = np.zeros((self.ndof, 1), dtype=MuscatFloat)
if K.nnz > 0 :
from Muscat.LinAlg.LinearSolver import LinearProblem
linSol = LinearProblem()
linSol.SetAlgo(self.linearSolver)
linSol.SetTolerance(self.tol)
linSol.SetOp(K)
linSol.u = self.u[self.free, 0]
self.u[self.free, 0] = linSol.Solve(rhs)
#else :
#print("'"+self.linearSolver + "' is not a valid linear solver")#pragma: no cover
#print('Please set a type of linear solver')#pragma: no cover
#raise Exception()#pragma: no cover
Debug('Post Process')
self.u = self.u + self.fixedValues
# Compute element elastic energy density
u_reshaped = self.u[self.edofMat]
nbBulkElements = self.mesh.GetNumberOfElements(self.mesh.GetElementsDimensionality())
u_reshaped.shape = (nbBulkElements, self.nodesPerElement*self.dofpernode)
Ku_reshaped = np.dot(u_reshaped, self.KE)
np.einsum('ij,ij->i', Ku_reshaped, u_reshaped, out=self.eed)
# we divide by the volume of one element to get the energy density
self.eed /= np.prod(self.mesh.props.get("spacing"))
Debug('Post Process Done')
[docs] def GenerateIJs(self):
# lazy generations of the IJ for the case when no density is given
Debug("GenerateIJs")
if self.iK is None:
nodesPerElement = 2**self.mesh.GetElementsDimensionality()
ones = np.ones((nodesPerElement*self.dofpernode, 1), dtype=MuscatIndex)
self.iK = np.kron(self.edofMat, ones).flatten()
ones.shape = (1, nodesPerElement*self.dofpernode)
self.jK = np.kron(self.edofMat, ones).flatten()
Debug("GenerateIJs Done")
[docs] def element_elastic_energy(self, Eeff= None, OnlyOnInterface = False):
if Eeff is None:
return self.eed
else:
nEeff = np.copy(Eeff)
bool_Eeff = Eeff<self.minthreshold
nEeff[bool_Eeff] = 0.0
return nEeff*self.eed
[docs] def nodal_elastic_energy(self, Eeff=None, OnlyOnInterface = False):
nbBulkElements = self.mesh.GetNumberOfElements(self.mesh.GetElementsDimensionality())
if Eeff is None:
Eeff = np.ones(nbBulkElements)
return node_averaged_element_field(self.element_elastic_energy(Eeff,OnlyOnInterface=OnlyOnInterface),self.mesh)
[docs] def Write(self):
if self.writer is not None:
self.writer.Write(self.mesh,
PointFields = [self.u, self.f],
PointFieldsNames= ["u", "f"]
)
[docs]def element_averaged_node_field(node_field,support):
nnodes = support.props.get("dimensions")
ndims = support.GetElementsDimensionality()
if ndims == 3:
result = np.zeros((nnodes[0]-1,nnodes[1]-1,nnodes[2]-1 ))
field = node_field.view()
field.shape = tuple(x for x in support.props.get("dimensions"))
result += field[0:-1, 0:-1,0:-1]
result += field[0:-1, 1: ,0:-1]
result += field[1: , 0:-1,0:-1]
result += field[1: , 1: ,0:-1]
result += field[0:-1, 0:-1,1: ]
result += field[0:-1, 1: ,1: ]
result += field[1: , 0:-1,1: ]
result += field[1: , 1: ,1: ]
return result/8
elif ndims == 2:
result = np.zeros((nnodes[0]-1,nnodes[1]-1))
field = node_field.view()
field.shape = tuple(x for x in support.props.get("dimensions"))
result += field[0:-1, 0:-1]
result += field[0:-1, 1: ]
result += field[1: , 0:-1]
result += field[1: , 1: ]
return result/4
else :
raise(Exception("only implemented for dim = 3 or 2"))
[docs]def node_averaged_element_field(element_field,support):
nnodes = support.props.get("dimensions")
ndims = support.GetElementsDimensionality()
if ndims == 3:
result = np.zeros((nnodes[0],nnodes[1],nnodes[2] ))
w = np.zeros((nnodes[0],nnodes[1],nnodes[2] ))
field = element_field.view()
field.shape = tuple(x-1 for x in nnodes)
result[0:-1, 0:-1,0:-1] += field
result[0:-1, 1: ,0:-1] += field
result[1: , 0:-1,0:-1] += field
result[1: , 1: ,0:-1] += field
result[0:-1, 0:-1,1: ] += field
result[0:-1, 1: ,1: ] += field
result[1: , 0:-1,1: ] += field
result[1: , 1: ,1: ] += field
w[0:-1, 0:-1,0:-1] += 1
w[0:-1, 1: ,0:-1] += 1
w[1: , 0:-1,0:-1] += 1
w[1: , 1: ,0:-1] += 1
w[0:-1, 0:-1,1: ] += 1
w[0:-1, 1: ,1: ] += 1
w[1: , 0:-1,1: ] += 1
w[1: , 1: ,1: ] += 1
return result/w
else:
result = np.zeros((nnodes[0],nnodes[1] ))
w = np.zeros((nnodes[0],nnodes[1] ))
field = element_field.view()
field.shape = tuple(x-1 for x in support.props.get("dimensions"))
result[0:-1, 0:-1] += field
result[0:-1, 1: ] += field
result[1: , 0:-1] += field
result[1: , 1: ] += field
w[0:-1, 0:-1] += 1
w[0:-1, 1: ] += 1
w[1: , 0:-1] += 1
w[1: , 1: ] += 1
return result/w
[docs]def CheckIntegrityThermal3D():
print('---------------------------- Thermal3D -------------------------------------------------')
from Muscat.Containers.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
import Muscat.IO.XdmfWriter as XdmfWriter
import time
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
nx = 11
ny = 12
nz = 13
myMesh = CreateConstantRectilinearMesh(dimensions=[nx,ny,nz], spacing=[0.1, 0.1, 10./(nz-1)], origin=[0,0,0])
dtag = myMesh.nodesTags.CreateTag("Block")
for x in range(nx):
for y in range(ny):
dtag.AddToTag(GetMonoIndex( [[x,y,0]],[nx,ny,nz])[0] )
print(myMesh)
# thermal problem
#dirichlet at plane z =0
myProblem = Fea()
myProblem.SetMesh(myMesh)
myProblem.AddDirichletCondition(NodeFilter(nTag="Block"), [True] , lambda x : [0] )
myProblem.AddNeumannCondition(NodeFilter(), [True] , lambda x : [1] )
starttime = time.time()
myElement = ElementaryMatrix()
myElement.physics = "thermal"
myElement.geoFactor = myMesh.props.get("spacing")
myProblem.BuildProblem(dofpernode = 1,
KOperator = myElement.GetTangentMatrix(),
MOperator = myElement.GetMassMatrix())
print("Time for Fea definition : " + str(time.time() -starttime))
myProblem.writer = XdmfWriter.XdmfWriter(TemporaryDirectory.GetTempPath()+"Test_constantRectilinearThermal.xdmf")
myProblem.writer.SetBinary()
myProblem.writer.SetTemporal()
myProblem.writer.Open()
myProblem.linearSolver = 'Direct'
myProblem.Solve()
myProblem.Write()
myProblem.writer.Close()
path = TemporaryDirectory.GetTempPath() +'TestThermal3D.xmf'
XdmfWriter.WriteMeshToXdmf(path,myMesh,
[myProblem.u, myProblem.f],
PointFieldsNames= ['Temperature','q'],
GridFieldsNames=[])
print('DONE')
print(TemporaryDirectory.GetTempPath())
print(np.max(myProblem.u))
if abs(np.max(myProblem.u)-50.) > 1e-4:
raise Exception(f'{abs(np.max(myProblem.u))} not equal to 50')# pragma: no cover
return("ok")
[docs]def CheckIntegrityDep3D():
import time
from Muscat.Containers.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
import Muscat.IO.XdmfWriter as XdmfWriter
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
print('------------------------------- Dep3D ----------------------------------------------')
nx = 11
ny = 12
nz = 13
myMesh = CreateConstantRectilinearMesh(dimensions=[nx,ny,nz],spacing=[0.1, 0.1, 10./(nz-1)], origin=[0, 0, 0] )
# block all the faces right
myProblem = Fea()
dirichletBcs = myProblem.dirichletBcs
dtagI = myMesh.nodesTags.CreateTag("BlockI")
dtagII = myMesh.nodesTags.CreateTag("BlockII")
for x in range(nx):
for y in range(ny):
dtagI.AddToTag(GetMonoIndex( [[x,y,0]],[nx,ny,nz])[0] )
dtagII.AddToTag(GetMonoIndex( [[x,y,nz-1]],[nx,ny,nz])[0] )
myProblem.AddDirichletCondition(NodeFilter(nTag="BlockI"), [True, True, True] , lambda x : [0,0,0] )
myProblem.AddDirichletCondition(NodeFilter(nTag="BlockII"), [True, True, True] , lambda x : [1,1,1] )
starttime = time.time()
myProblem.SetMesh(myMesh)
myProblem.BuildProblem(dofpernode = 3)
print("Time for Fea definition : " + str(time.time() -starttime))
starttime = time.time()
densities = np.ones(myMesh.GetNumberOfElements(dim=3))
myProblem.Solve(densities)
print("Time for Fea solve : " + str(time.time() -starttime))
XdmfWriter.WriteMeshToXdmf(TemporaryDirectory.GetTempPath() +'TestDep3D.xmf',myMesh,
[myProblem.u, myProblem.f, myProblem.nodal_elastic_energy() ],
[densities, myProblem.element_elastic_energy() ] ,
[],
PointFieldsNames= ['Dep','F','NEnergie'],
CellFieldsNames=['densities','EEnergie'],
GridFieldsNames=[])
print(np.max(myProblem.u))
if abs(np.max(myProblem.u)-1.00215295) > 1e-5:
print(TemporaryDirectory.GetTempPath())
raise Exception("The value must be 1.00215295")# pragma: no cover
[docs]def CheckIntegrityThermal2D():
from Muscat.Containers.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
import Muscat.IO.XdmfWriter as XdmfWriter
import time
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
print('---------------------------- Thermal2D -------------------------------------------------')
nx = 11
ny = 11
myMesh = CreateConstantRectilinearMesh(dimensions=[nx,ny],spacing=[0.1, 0.1], origin=[0, 0] )
dtagI = myMesh.nodesTags.CreateTag("Block")
dtagII = myMesh.nodesTags.CreateTag("Load")
for x in range(nx):
dtagI.AddToTag(GetMonoIndex( [[x,0]],[nx,ny])[0] )
for y in range(ny):
dtagII.AddToTag(GetMonoIndex([[x,y]],[nx,ny])[0] )
print(myMesh)
myProblem = Fea()
myProblem.SetMesh(myMesh)
myProblem.AddDirichletCondition(NodeFilter(nTag="Block"), [True],lambda x : [0] )
myProblem.AddNeumannCondition(NodeFilter(nTag="Load"), [True], lambda x : [1] )
starttime = time.time()
myElement = ElementaryMatrix(dim=2,physics="thermal")
myElement.geoFactor = myMesh.props.get("spacing")
myProblem.BuildProblem(dofpernode = 1,
KOperator = myElement.GetTangentMatrix(),
MOperator = myElement.GetMassMatrix())
print("Time for Fea definition : " + str(time.time() -starttime))
myProblem.writer = XdmfWriter.XdmfWriter(TemporaryDirectory.GetTempPath()+"TestProblemWriterThermal2D.xdmf")
myProblem.writer.SetBinary()
myProblem.writer.SetTemporal()
myProblem.writer.Open()
myProblem.linearSolver = 'Direct'
myProblem.Solve()
myProblem.Write()
myProblem.element_elastic_energy()
myProblem.Write()
myProblem.writer.Close()
path = TemporaryDirectory.GetTempPath() +'TestThermal2D.xmf'
XdmfWriter.WriteMeshToXdmf(path,myMesh,
[myProblem.u, myProblem.f],
PointFieldsNames= ['Temperature','q'],
GridFieldsNames=[])
print('DONE')
print(np.max(myProblem.u))
if abs(np.max(myProblem.u)-.5) > 1e-5:
raise Exception(f"{np.max(myProblem.u)} grater than 5")# pragma: no cover
return 'OK'
[docs]def CheckIntegrityDep2D():
from Muscat.Containers.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
import Muscat.IO.XdmfWriter as XdmfWriter
import time
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
print('----------------------- 2D dep ------------------------------------------------------')
nx = 11
ny = 12
myMesh = CreateConstantRectilinearMesh(dimensions=[nx,ny], spacing=[1./(nx-1), 1./(ny-1)], origin=[0, 0] )
dtagI = myMesh.nodesTags.CreateTag("BlockI")
dtagII = myMesh.nodesTags.CreateTag("BlockII")
for x in range(nx):
dtagI.AddToTag(GetMonoIndex( [[x,0]],[nx,ny])[0] )
dtagII.AddToTag(GetMonoIndex([[x,ny-1]],[nx,ny])[0] )
print(myMesh)
myProblem = Fea()
myProblem.SetMesh(myMesh)
myProblem.AddDirichletCondition(NodeFilter(nTag="BlockI"), [True, True], lambda x : [0, 0] )
myProblem.AddDirichletCondition(NodeFilter(nTag="BlockII"), [True, True], lambda x : [1, 1] )
starttime = time.time()
myProblem.BuildProblem(dofpernode = 2)
print("Time for Fea definition : " + str(time.time() -starttime))
starttime = time.time()
densities = np.ones(myMesh.GetNumberOfElements(dim=2))
myProblem.Solve(densities)
print("Time for Fea solve : " + str(time.time() -starttime))
print(myProblem.u.T)
# build mass matrix
myProblem.BuildMassMatrix()
myProblem.BuildMassMatrix(densities)
XdmfWriter.WriteMeshToXdmf(TemporaryDirectory.GetTempPath() +'TestDep2D.xmf',myMesh,
[myProblem.u, myProblem.f, myProblem.nodal_elastic_energy(), myProblem.fixed.astype(int) ],
[densities, myProblem.element_elastic_energy() ] ,
[],
PointFieldsNames= ['Dep','F','NEnergie', 'ufixed'],
CellFieldsNames=['densities','EEnergie'],
GridFieldsNames=[])
return 'ok'
[docs]def CheckIntegrity():
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
print(CheckIntegrityThermal3D())
print(CheckIntegrityDep3D())
print(CheckIntegrityThermal2D())
print(CheckIntegrityDep2D())
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
print(TemporaryDirectory.GetTempPath())
return "ok"
if __name__ == '__main__':
print(CheckIntegrity()) #pragma: no cover