#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Dec 13 15:10:47 2023
@author: nchevaug
"""
from typing import cast
import numpy as np
import scipy as sp
try:
import sksparse.cholmod
except :
have_cholmod = False
from pyxfem.core.config import npfloat
import pyxfem.core.function_spaces as fs
import pyxfem.core.mapping as mapping
import pyxfem.core.quadrature as quad
import pyxfem.core.sparse_tools as st
import pyxfem.core.level_set_mesh_cut as lscut
#import pyxfem.xmesh
import pyxfem.core.gmsh_mesh_model_interface as mi
import pyxfem.core.mesh_simplex_one_level as sm
import pyxfem.core.sparse_tools as st
[docs]
def getDofsAndInterpolationPoints(modelName, groupeName2DimPhyTag, groupeNames, space) :
dirichlet2dofs = {}
dirichlet2xy = {}
for physName in groupeNames :
dim, phystag = groupeName2DimPhyTag[physName]
if space.interpolatoryAt == 'vertex':
nodes_tag, coord = mi.get_vertices_on_phys_tag(modelName, dim, phystag)
dofsid = space.vertexid2dofid(nodes_tag)
dirichlet2dofs[physName] = dofsid
dirichlet2xy[physName] = coord.reshape((-1,2))
elif space.interpolatoryAt == 'midEdge' :
e2v_phys = mi.get_toplevel_elements_vertices_on_phys_tag(modelName, dim, phystag).get('edge', np.zeros((0,2), dtype = np.int64 ))
eid_space, exy = space.edge2edgeid(e2v_phys, getxy = True)
dirichlet2dofs[physName] = space.edgeid2dofid(eid_space).squeeze()
dirichlet2xy[physName] = 0.5*np.sum(exy, axis = 1)
else :
raise ValueError('Error : ', getDofsAndInterpolationPoints.__name__, ' not coded for space.interpolatoryAt ', space.interpolatoryAt, '.' )
return dirichlet2dofs, dirichlet2xy
[docs]
class ThermicSolver():
[docs]
@staticmethod
def setconductivity_fromls(conductivities, cutmesh, ls_tri, ls_tri_cut, rmap):
if np.isscalar(conductivities) : conductivity_in, conductivity_out = [conductivities]*2
else : conductivity_in, conductivity_out = conductivities
if cutmesh is None : return conductivities
conductivity_eval_tri = np.zeros(rmap.nbevals)
side_appro = lscut.compute_face_side(ls_tri)
side_eval = lscut.compute_face_side(ls_tri_cut)
conductivity_eval_tri[rmap.noson2eval_index] = np.where( side_appro[rmap.approEQeval] == -1, conductivity_in, conductivity_out)
conductivity_eval_tri[rmap.son2eval_index] = np.where( side_eval[rmap.sons] == -1, conductivity_in, conductivity_out)
return conductivity_eval_tri
[docs]
@staticmethod
def setconductivity_fromls0(conductivities, ls_tri):
if np.isscalar(conductivities) : conductivity_in, conductivity_out = [conductivities]*2
else : conductivity_in, conductivity_out = conductivities
side = lscut.compute_face_side(ls_tri)
conductivity = np.where( side == -1, conductivity_in, conductivity_out)
return conductivity
def __init__(self, modelName, physical_dict, conductivity, dirichlet, source = 0., space = 'P1C', lsfunforconductivity = None, enrich = None, xmesh_relay = False, dosvd = False):
''' Thermic problem constructor
input :
- modelName : a string corresponding to the name of the model to look at in gmsh.
Used for getting the mesh and informations on the classification of the mesh entities
- physical_dict : a dictionnary, giving for each named physcical group the associated (dim, tag) pair defining the associated Group in the gmsh model,
where dim is the dimension dimension of the physical entity and tag its tag.
- conductivity : define the conductivity. either globally or relative to the leveset values lsfun for conductivity
- dirichlet : a map associating the name of a physical group to a Dirichlet value or function of xy
- source : the source term value or a function of xy to compute it at point
- space : base finite element space name. For Now, only 'P1C' (linear continuous element) and 'P1NC' (linear discontinuous Crouzeix-RAviart Element)
- lsfunforconductivity: None or a function of space whose sign decide the material conductivity ( ls < 0 -> k <- conductivity[0], ls >=0 -> k <- conductivity[1]
- enrich : for now, one of None, 'xRidgeSpaceP1', 'xFEMspace'. if not None, the space will be enriched to allow for gradient discontinuity alng iso-0 of lsfunfor conductivity
- xmesh_relay: False, or True : an ungly flagto set if we do xmesh relaying algorithm or not to be conform to the interface defined by lsfunforconductivity
- dosvd : False : an other ugly tag to compute or not an svd of the finiteelement space .... This was used for debugging purpose when testing the new xRidgeSpaceP1 enrichment
'''
if enrich and lsfunforconductivity is None :
raise ValueError("can't enrich if no levelset function is given")
self.modelName = modelName
self.physical_dict = physical_dict
self.dirichlet = dirichlet
# if xmesh_relay :
# raise NotImplementedError("We need to reimplement xmesh for this to work")
# m = xmesh.Mesh(modelName)
# xy0 = m.xy
# topo = xmesh.SimplexTopo(m.tri2node)
# tris = topo.tri2node
# self.xy = xy0
# self.tris = tris
# if lsfunforconductivity :
# phi = lsfunforconductivity
# ls = phi(self.xy)
# xy, front, sol = xmesh.initialize_front(topo, xy0, phi)
# self.xy = xy
# else:
self.xy = mi.get_vertices_coord(modelName)
self.tris, __ = mi.get_triangles(modelName)
integrationRule0 = quad.T3_gauss(0)
self.rmap = None
if lsfunforconductivity :
phi = lsfunforconductivity
ls = phi(self.xy)
lscut.fit_to_vertices(ls, self.tris, absfittol = 1.e-12, relfittol = 1.e-3)
if xmesh_relay :
#nf = len(self.tris)
#hmin = np.sqrt(4/nf)
self.tris_xy = self.xy[self.tris]
self.tris_map = mapping.MappingT3(self.tris_xy, hmin =0., alpha = 0.01, p = 1)
self.space = fs.FEMSpaceP1("Temperature", sm.SMesh(self.xy, self.tris), vecdim = 1, tris_map = self.tris_map )
self.B0 = self.space.operator_dof2grad_tri(integrationRule0.uv, limited = False)
self.B = self.space.operator_dof2grad_tri(integrationRule0.uv, limited = True)
#print('diff B', np.linalg.norm(self.B0.data - self.B.data))
self.N = self.space.operator_dof2val_tri(quad.T3_gauss(0).uv)
J = self.tris_map.J
self.J = J
conductivities = ThermicSolver.setconductivity_fromls0(conductivity, ls[self.tris])
D = (J[:, np.newaxis]* np.tensordot(conductivities, np.eye(2).flatten()*integrationRule0.w[0], axes = 0 )).reshape((-1, 2, 2))
self.D = st.regular_block_diag_to_csr(D)
self.NatVertices = self.space.operator_dof2val_tri(quad.T3_nodes().uv)
self.Fint = integrationRule0.w[0]*source *self.N.T.dot(J)
else :
cut_mesh, cut_ls = lscut.levelset_isocut(self.xy, self.tris, ls, returnparent = True)
rmap = fs.meshcutdata2remap_helper(cut_mesh, hmin = 0., alpha = 0.1, p = 1)
self.tris_xy = rmap.get_eval_tris()
self.space = fs.FEMSpaceP1("Temperature", sm.SMesh(self.xy, self.tris), vecdim = 1, tris_map = rmap.parent_mapping)
self.B = self.space.operator_dof2grad_tri(integrationRule0.uv,rmap)
self.B0 = self.B
self.N = self.space.operator_dof2val_tri(quad.T3_gauss(0).uv, rmap)
self.NatVertices = self.space.operator_dof2val_tri(quad.T3_nodes().uv, rmap)
if enrich :
spaceP1 = fs.FEMSpaceP1("Temperature", sm.SMesh(self.xy, self.tris), vecdim = 1, tris_map = rmap.parent_mapping)
if enrich == 'xFEMSpace' :
self.xspace = fs.XFEMSpace('ridge', fs.FEMSpaceP1, fs.Ridge, cut_mesh, ls, cut_ls, rmap = rmap)
elif enrich == 'xRidgeSpaceP1':
self.xspace = fs.XRidgeSpaceP1('ridge',cut_mesh, ls, cut_ls, rmap = rmap, dosvd = False)
else :
raise ValueError(f"enrichement {enrich} is not known")
self.space = fs.SumSpace([spaceP1, self.xspace])
self.B = self.space.operator_dof2grad_tri(integrationRule0.uv, rmap)
self.B0 = self.B
self.N = self.space.operator_dof2val_tri(quad.T3_gauss(0).uv, rmap)
self.NatVertices = self.space.operator_dof2val_tri(quad.T3_nodes().uv, rmap)
neval = rmap.nbevals
J = np.zeros(neval)
J[rmap.noson2eval_index] = rmap.parent_mapping.J[rmap.approEQeval]
J[rmap.son2eval_index] = rmap.son_mapping.J[rmap.sons]
self.J = J
#print(np.min(J), np.max(J))
conductivities = ThermicSolver.setconductivity_fromls(conductivity, cut_mesh, ls_tri=ls[self.tris], ls_tri_cut= cut_ls[cut_mesh.new_tris], rmap= rmap)
conductivities = np.asarray(conductivities, npfloat)
ws = np.asarray(np.eye(2).flatten()*integrationRule0.w[0], npfloat)
tmp = np.tensordot(conductivities, ws, axes = 0 )
Dval = (J[:, np.newaxis]*tmp).flatten()
Dj = np.repeat(np.arange(2*neval).reshape((-1,2)), 2, axis = 0).flatten()
Di = np.arange(2*neval).repeat(2)
self.D = sp.sparse.csr_array((Dval, (Di, Dj)), (2*neval, 2*neval) )
self.Fint = integrationRule0.w[0]*source *self.N.T.dot(J)
self.rmap = rmap
else :
self.tris_xy = self.xy[self.tris]
self.tris_map = mapping.MappingT3(self.tris_xy)
if space =='P1C':
self.space = fs.FEMSpaceP1("Temperature", sm.SMesh(self.xy, self.tris), vecdim = 1, tris_map = self.tris_map)
elif space =='P1NC':
self.space = fs.FEMSpaceP1NC("Temperature", sm.SMesh(self.xy, self.tris), vecdim = 1, tris_map = self.tris_map)
else :
print('Error : space ', space, ' not coded in ', ThermicSolver.__name__,'.' )
raise NotImplementedError(f"space not coded")
self.B = self.space.operator_dof2grad_tri(integrationRule0.uv)
self.B0 = self.B
nt = len(self.tris)
J = self.tris_map.J
self.J = J
Dval = np.tensordot(J, np.eye(2).flatten()*conductivity*integrationRule0.w[0], axes = 0).flatten()
Dj = np.repeat(np.arange(2*nt).reshape((-1,2)), 2, axis = 0).flatten()
Di = np.arange(2*nt).repeat(2)
self.D = sp.sparse.csr_array((Dval, (Di, Dj)), (2*nt, 2*nt) )
self.N = self.space.operator_dof2val_tri(integrationRule0.uv)
self.NatVertices = self.space.operator_dof2val_tri(quad.T3_nodes().uv)
if callable(source) :
quadrature = quad.T3_gauss(5)
npt = len(quadrature.w)
uv = quadrature.uv
weight = quadrature.w
NP1 = fs.FEMSpaceP1.N(uv)
tris_gp_xy = np.einsum('ijk, ikl-> ijl', NP1[np.newaxis, :,:], self.tris_xy)
sourcegp = source(tris_gp_xy)
wsourcegp = weight[np.newaxis,:] *J[:, np.newaxis]* sourcegp
Nsource = self.space.operator_dof2val_tri(uv)
self.Fint = Nsource.T.dot(wsourcegp.flatten())
elif np.isscalar(source):
self.Fint = integrationRule0.w[0]*source *self.N.T.dot(J)
else :
raise ValueError("the source term must be a scalar or be evaluatable")
[docs]
def solve(self, solver = "default", compute_rcond = False):
if type(self.space) == fs.SumSpace : space = self.space.list_of_space[0]
else : space = self.space
self.dirichlet2dofs, self.dirichlet2xy = getDofsAndInterpolationPoints(self.modelName, self.physical_dict, self.dirichlet.keys(), space)
fixeddofs = np.hstack([ np.array(dofs, dtype = np.int64) for dofs in self.dirichlet2dofs.values() ])
fixedvalues = []
for xy, dofs, dval in zip(self.dirichlet2xy.values(), self.dirichlet2dofs.values(), self.dirichlet.values()):
if np.isscalar(dval) :
dval = cast(npfloat, dval)
fixedvalues.append(np.ones(len(dofs))*dval)
elif callable(dval) :
fixedvalues.append(dval(xy))
else :
raise ValueError("val mus be a scalar or a callable")
fixedvalues = np.hstack(fixedvalues)
freedof = np.setdiff1d(np.arange(self.B.shape[1]), fixeddofs)
Bfree = self.B[:, freedof]
U = np.zeros(self.B.shape[1])
U[fixeddofs] = fixedvalues
gradTdiri = self.B[:, fixeddofs].dot(fixedvalues)
Fdiri = - Bfree.T.dot(self.D.dot(gradTdiri))
F = self.Fint[freedof] + Fdiri
if solver == "default" :
Kfree = Bfree.T.dot(self.D).dot(Bfree)
luKfree = sp.sparse.linalg.splu(Kfree)
U[freedof] = luKfree.solve(F)
if compute_rcond :
rcond = st.rcond_estimate(Kfree, luKfree)
return U, rcond
return U
if solver == "cholmod":
Kfree = Bfree.T.dot(self.D).dot(Bfree)
chol = sksparse.cholmod.cholesky(Kfree)
U[freedof] = chol(F)
if compute_rcond:
diag = chol.D()
rcond = (np.min(diag)/np.max(diag))**2
return U, rcond
return U
if solver == "AATcholmod":
sqrtD = sp.sparse.csr_array( (np.sqrt(self.D.data), self.D.indices, self.D.indptr), self.D.shape)
chol = sksparse.cholmod.cholesky_AAt( (sqrtD.dot(Bfree)).T)
U[freedof] = chol(F)
if compute_rcond :
raise NotImplementedError('rcond is not implemeted for AATcholmod')
return U
print('solver '+ solver + ' is not known !')
raise ValueError('solver can only be one of "default", "cholmod", "AATcholmod"' )
[docs]
def postproVal(self, U):
return self.tris_xy, self.NatVertices.dot(U).reshape((-1, 3))
[docs]
def postproGrad(self, U):
return self.tris_xy, self.B.dot(U).reshape((-1, 2))
[docs]
def H1_error(self, U, gradexactfun, quadrature = quad.T3_gauss(0)):
npt = len(quadrature.w)
uv = quadrature.uv
weight = quadrature.w
B = self.space.operator_dof2grad_tri(uv, self.rmap)
gradNum = B.dot(U).reshape((-1, npt, 2))
NP1 = fs.FEMSpaceP1.N(uv)
tris_gp_xy = np.einsum('ijk, ikl-> ijl', NP1[np.newaxis, :,:], self.tris_xy)
gradExact = gradexactfun(tris_gp_xy)
error2_gp = np.linalg.norm((gradNum - gradExact), axis = 2)**2
error2 = np.sum(self.J*np.einsum( 'ij, ij-> i', weight[np.newaxis, :] , error2_gp))
return np.sqrt(error2)
[docs]
def maxValError(self, U, exactfun, samplingpoints = quad.T3_gauss(0).uv) :
uv = samplingpoints
npt = len(uv)
N = self.space.operator_dof2val_tri(uv, self.rmap)
valNum = N.dot(U).reshape((-1, npt))
NP1 = fs.FEMSpaceP1.N(uv)
tris_gp_xy = np.einsum('ijk, ikl-> ijl', NP1[np.newaxis, :,:], self.tris_xy)
valExact = exactfun(tris_gp_xy)
error_gp = np.abs(valNum - valExact)
return np.max(error_gp)
[docs]
def maxGradError(self, U, gradexactfun, samplingpoints = quad.T3_gauss(0).uv) :
uv = samplingpoints
npt = len(uv)
B = self.space.operator_dof2grad_tri(uv, self.rmap)
gradNum = B.dot(U).reshape((-1, npt, 2))
NP1 = fs.FEMSpaceP1.N(uv)
tris_gp_xy = np.einsum('ijk, ikl-> ijl', NP1[np.newaxis, :,:], self.tris_xy)
gradExact = gradexactfun(tris_gp_xy)
error_gp = np.linalg.norm((gradNum - gradExact), axis = 2)
return np.max(error_gp)
[docs]
def L2_error(self, U, exactfun, quadrature = quad.T3_gauss(0)):
npt = len(quadrature.w)
uv = quadrature.uv
weight = quadrature.w
N = self.space.operator_dof2val_tri(uv, self.rmap)
ValNum = N.dot(U).reshape((-1, npt))
NP1 = fs.FEMSpaceP1.N(uv)
tris_gp_xy = np.einsum('ijk, ikl-> ijl', NP1[np.newaxis, :,:], self.tris_xy)
ValExact = exactfun(tris_gp_xy)
error2_gp = (ValNum - ValExact)**2
error2 = np.sum(self.J*np.einsum( 'ij, ij-> i', weight[np.newaxis, :] , error2_gp))
return np.sqrt(error2)
if __name__== '__main__':
import gmsh
from pyxfem.core.exec_tools import Profile
import pyxfem.core.gmsh_mesh_generator as mg
import pyxfem.postpro.gmsh_post_pro_interface as gp
elementSize = 0.2
solver = "default" #"cholmod" # "default", "AATcholmod"
popUpGmsh = True
Profile.doprofile = False
dofem = ['P1C', 'P1NC']
dofem2conductivity = True
doxfem2conductivity = ['xFEMSpace', 'xRidgeSpaceP1']
doxfem_bench = True
# dofem = ['P1C']
# dofem2conductivity = False
# doxfem2conductivity = []
# doxfem_bench = False
gmsh.initialize()
gmsh.option.setNumber('PostProcessing.AnimationCycle', 1)
gmsh.option.setNumber('Mesh.SurfaceEdges', 1 if elementSize > 0.1 else 0)
modelName, physical_dict = mg.l_shape(elementSize)
nv = mi.nb_vertices(modelName)
nf = mi.nb_triangles(modelName)
print('nb elements =', nf)
print('nb vertices =', nv)
physname_to_dirichlet_value = {'dOmega_left':0.33, 'dOmega_right' :0.1}
conductivity = 1.
for spaceName in dofem :
print('Solving a fem problem using space ', spaceName )
pbfem = ThermicSolver(modelName, physical_dict, conductivity, physname_to_dirichlet_value, space = spaceName, enrich = False )
U = pbfem.solve(solver = solver)
trisxy, Ttris = pbfem.postproVal(U)
trisxy, GTtris = pbfem.postproGrad(U)
gp.listPlotFieldTri(trisxy, Ttris, P0 = False, viewname ='T_'+ spaceName, NbIso = 20)
gp.listPlotFieldTri(trisxy, GTtris, P0 = True, viewname ='GradT_'+ spaceName, NbIso = 20, VectorType = "Arrow")
space = pbfem.space
if type(space) == fs.FEMSpaceP1NC:
N_edge = space.operator_edgetrace(quad.Edge_nodes())
valEdge = N_edge.dot(U).reshape(-1,2)
gp.listPlotFieldLine(space.m.e2v_coord(), valEdge[:,0], P0 = False, viewname ='T_EdgeL', LineType = "3D cylinder", LineWidth = 4)
gp.listPlotFieldLine(space.m.e2v_coord(), valEdge[:,1], P0 = False, viewname ='T_EdgeR', LineType = "3D cylinder", LineWidth = 4)
if dofem2conductivity :
print('Solving a fem problem with 2 conductivities defined by a levelset')
lsfun = lscut.SimpleLevelSetFunction.disk([-0.7,0.5], 1.)
pbfem = ThermicSolver(modelName, physical_dict, [20*conductivity, conductivity], physname_to_dirichlet_value, lsfunforconductivity = lsfun, enrich = False )
U = pbfem.solve(solver = solver)
trisxy, Ttris = pbfem.postproVal(U)
trisxy, GTtris = pbfem.postproGrad(U)
gp.listPlotFieldTri(trisxy, Ttris, P0 = False, viewname ='T2cond_fem', NbIso = 20)
gp.listPlotFieldTri(trisxy, GTtris, P0 = True, viewname ='GradT2cond_fem', NbIso = 20, VectorType = "Displacement")
for enrichmentName in doxfem2conductivity :
print('Solving a xfem problem with 2 conductivities defined by a levelset. enrichment is ', enrichmentName)
lsfun = lscut.SimpleLevelSetFunction.disk([-0.7,0.5], 1.)
pbxfem = ThermicSolver(modelName, physical_dict, [20*conductivity, conductivity], physname_to_dirichlet_value, lsfunforconductivity = lsfun, enrich = enrichmentName)
print(" nb enriched dofs = ", pbxfem.xspace.scalar_size())
U = pbxfem.solve(solver = solver)
trisxy, Ttris = pbxfem.postproVal(U)
trisxy, GTtris = pbxfem.postproGrad(U)
gp.listPlotFieldTri(trisxy, Ttris, P0 = False, viewname ='T2cond_'+enrichmentName, NbIso = 20)
gp.listPlotFieldTri(trisxy, GTtris, P0 = True, viewname ='GradT2cond_'+enrichmentName, NbIso = 20, VectorType = "Displacement")
if doxfem_bench :
print('Solving a xfem problem with 2 conductivities defined by a more complicated levelset.')
x = np.linspace(-0.8, 0.8, 4)
y = np.linspace(-0.8, 0.8, 4)
lsfun = lscut.SimpleLevelSetFunction.union([lscut.SimpleLevelSetFunction.disk([xi,yi], 0.10) for xi in x for yi in y] )
pbxfem = ThermicSolver(modelName, physical_dict, [conductivity*0.001, conductivity], physname_to_dirichlet_value, lsfunforconductivity = lsfun, enrich = 'xFEMSpace')
print(" nb enriched dofs = ", pbxfem.xspace.scalar_size())
U = pbxfem.solve(solver = solver)
trisxy, Ttris = pbxfem.postproVal(U)
trisxy, GTtris = pbxfem.postproGrad(U)
gp.listPlotFieldTri(trisxy, Ttris, P0 = False, viewname ='BenchXfem_T', NbIso = 20)
gp.listPlotFieldTri(trisxy, GTtris, P0 = True, viewname ='BenchXfem_GradT', NbIso = 20, VectorType = "Displacement")
if popUpGmsh : gmsh.fltk.run()
gmsh.finalize()
Profile.print_stats()