Coverage for pyxfem/thermic/thermic_solver.py: 0%
314 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3"""
4Created on Wed Dec 13 15:10:47 2023
6@author: nchevaug
7"""
8from typing import cast
9import numpy as np
10import scipy as sp
11try:
12 import sksparse.cholmod
13except :
14 have_cholmod = False
16from pyxfem.core.config import npfloat
17import pyxfem.core.function_spaces as fs
18import pyxfem.core.mapping as mapping
19import pyxfem.core.quadrature as quad
20import pyxfem.core.sparse_tools as st
21import pyxfem.core.level_set_mesh_cut as lscut
22#import pyxfem.xmesh
23import pyxfem.core.gmsh_mesh_model_interface as mi
24import pyxfem.core.mesh_simplex_one_level as sm
25import pyxfem.core.sparse_tools as st
27def getDofsAndInterpolationPoints(modelName, groupeName2DimPhyTag, groupeNames, space) :
28 dirichlet2dofs = {}
29 dirichlet2xy = {}
30 for physName in groupeNames :
31 dim, phystag = groupeName2DimPhyTag[physName]
32 if space.interpolatoryAt == 'vertex':
33 nodes_tag, coord = mi.get_vertices_on_phys_tag(modelName, dim, phystag)
34 dofsid = space.vertexid2dofid(nodes_tag)
35 dirichlet2dofs[physName] = dofsid
36 dirichlet2xy[physName] = coord.reshape((-1,2))
37 elif space.interpolatoryAt == 'midEdge' :
38 e2v_phys = mi.get_toplevel_elements_vertices_on_phys_tag(modelName, dim, phystag).get('edge', np.zeros((0,2), dtype = np.int64 ))
39 eid_space, exy = space.edge2edgeid(e2v_phys, getxy = True)
40 dirichlet2dofs[physName] = space.edgeid2dofid(eid_space).squeeze()
41 dirichlet2xy[physName] = 0.5*np.sum(exy, axis = 1)
42 else :
43 raise ValueError('Error : ', getDofsAndInterpolationPoints.__name__, ' not coded for space.interpolatoryAt ', space.interpolatoryAt, '.' )
44 return dirichlet2dofs, dirichlet2xy
46class ThermicSolver():
47 @staticmethod
48 def setconductivity_fromls(conductivities, cutmesh, ls_tri, ls_tri_cut, rmap):
49 if np.isscalar(conductivities) : conductivity_in, conductivity_out = [conductivities]*2
50 else : conductivity_in, conductivity_out = conductivities
51 if cutmesh is None : return conductivities
52 conductivity_eval_tri = np.zeros(rmap.nbevals)
53 side_appro = lscut.compute_face_side(ls_tri)
54 side_eval = lscut.compute_face_side(ls_tri_cut)
55 conductivity_eval_tri[rmap.noson2eval_index] = np.where( side_appro[rmap.approEQeval] == -1, conductivity_in, conductivity_out)
56 conductivity_eval_tri[rmap.son2eval_index] = np.where( side_eval[rmap.sons] == -1, conductivity_in, conductivity_out)
57 return conductivity_eval_tri
58 @staticmethod
59 def setconductivity_fromls0(conductivities, ls_tri):
60 if np.isscalar(conductivities) : conductivity_in, conductivity_out = [conductivities]*2
61 else : conductivity_in, conductivity_out = conductivities
62 side = lscut.compute_face_side(ls_tri)
63 conductivity = np.where( side == -1, conductivity_in, conductivity_out)
64 return conductivity
66 def __init__(self, modelName, physical_dict, conductivity, dirichlet, source = 0., space = 'P1C', lsfunforconductivity = None, enrich = None, xmesh_relay = False, dosvd = False):
67 ''' Thermic problem constructor
68 input :
69 - modelName : a string corresponding to the name of the model to look at in gmsh.
70 Used for getting the mesh and informations on the classification of the mesh entities
71 - physical_dict : a dictionnary, giving for each named physcical group the associated (dim, tag) pair defining the associated Group in the gmsh model,
72 where dim is the dimension dimension of the physical entity and tag its tag.
73 - conductivity : define the conductivity. either globally or relative to the leveset values lsfun for conductivity
74 - dirichlet : a map associating the name of a physical group to a Dirichlet value or function of xy
75 - source : the source term value or a function of xy to compute it at point
76 - space : base finite element space name. For Now, only 'P1C' (linear continuous element) and 'P1NC' (linear discontinuous Crouzeix-RAviart Element)
77 - lsfunforconductivity: None or a function of space whose sign decide the material conductivity ( ls < 0 -> k <- conductivity[0], ls >=0 -> k <- conductivity[1]
78 - 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
79 - 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
80 - 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
81 '''
82 if enrich and lsfunforconductivity is None :
83 raise ValueError("can't enrich if no levelset function is given")
84 self.modelName = modelName
85 self.physical_dict = physical_dict
86 self.dirichlet = dirichlet
87 # if xmesh_relay :
88 # raise NotImplementedError("We need to reimplement xmesh for this to work")
89 # m = xmesh.Mesh(modelName)
90 # xy0 = m.xy
91 # topo = xmesh.SimplexTopo(m.tri2node)
92 # tris = topo.tri2node
93 # self.xy = xy0
94 # self.tris = tris
95 # if lsfunforconductivity :
96 # phi = lsfunforconductivity
97 # ls = phi(self.xy)
98 # xy, front, sol = xmesh.initialize_front(topo, xy0, phi)
99 # self.xy = xy
100 # else:
101 self.xy = mi.get_vertices_coord(modelName)
102 self.tris, __ = mi.get_triangles(modelName)
103 integrationRule0 = quad.T3_gauss(0)
104 self.rmap = None
105 if lsfunforconductivity :
106 phi = lsfunforconductivity
107 ls = phi(self.xy)
108 lscut.fit_to_vertices(ls, self.tris, absfittol = 1.e-12, relfittol = 1.e-3)
109 if xmesh_relay :
110 #nf = len(self.tris)
111 #hmin = np.sqrt(4/nf)
112 self.tris_xy = self.xy[self.tris]
113 self.tris_map = mapping.MappingT3(self.tris_xy, hmin =0., alpha = 0.01, p = 1)
114 self.space = fs.FEMSpaceP1("Temperature", sm.SMesh(self.xy, self.tris), vecdim = 1, tris_map = self.tris_map )
115 self.B0 = self.space.operator_dof2grad_tri(integrationRule0.uv, limited = False)
116 self.B = self.space.operator_dof2grad_tri(integrationRule0.uv, limited = True)
117 #print('diff B', np.linalg.norm(self.B0.data - self.B.data))
118 self.N = self.space.operator_dof2val_tri(quad.T3_gauss(0).uv)
119 J = self.tris_map.J
120 self.J = J
121 conductivities = ThermicSolver.setconductivity_fromls0(conductivity, ls[self.tris])
122 D = (J[:, np.newaxis]* np.tensordot(conductivities, np.eye(2).flatten()*integrationRule0.w[0], axes = 0 )).reshape((-1, 2, 2))
123 self.D = st.regular_block_diag_to_csr(D)
124 self.NatVertices = self.space.operator_dof2val_tri(quad.T3_nodes().uv)
126 self.Fint = integrationRule0.w[0]*source *self.N.T.dot(J)
127 else :
128 cut_mesh, cut_ls = lscut.levelset_isocut(self.xy, self.tris, ls, returnparent = True)
129 rmap = fs.meshcutdata2remap_helper(cut_mesh, hmin = 0., alpha = 0.1, p = 1)
130 self.tris_xy = rmap.get_eval_tris()
131 self.space = fs.FEMSpaceP1("Temperature", sm.SMesh(self.xy, self.tris), vecdim = 1, tris_map = rmap.parent_mapping)
132 self.B = self.space.operator_dof2grad_tri(integrationRule0.uv,rmap)
133 self.B0 = self.B
134 self.N = self.space.operator_dof2val_tri(quad.T3_gauss(0).uv, rmap)
135 self.NatVertices = self.space.operator_dof2val_tri(quad.T3_nodes().uv, rmap)
137 if enrich :
138 spaceP1 = fs.FEMSpaceP1("Temperature", sm.SMesh(self.xy, self.tris), vecdim = 1, tris_map = rmap.parent_mapping)
139 if enrich == 'xFEMSpace' :
140 self.xspace = fs.XFEMSpace('ridge', fs.FEMSpaceP1, fs.Ridge, cut_mesh, ls, cut_ls, rmap = rmap)
141 elif enrich == 'xRidgeSpaceP1':
142 self.xspace = fs.XRidgeSpaceP1('ridge',cut_mesh, ls, cut_ls, rmap = rmap, dosvd = False)
143 else :
144 raise ValueError(f"enrichement {enrich} is not known")
145 self.space = fs.SumSpace([spaceP1, self.xspace])
146 self.B = self.space.operator_dof2grad_tri(integrationRule0.uv, rmap)
147 self.B0 = self.B
148 self.N = self.space.operator_dof2val_tri(quad.T3_gauss(0).uv, rmap)
149 self.NatVertices = self.space.operator_dof2val_tri(quad.T3_nodes().uv, rmap)
150 neval = rmap.nbevals
151 J = np.zeros(neval)
152 J[rmap.noson2eval_index] = rmap.parent_mapping.J[rmap.approEQeval]
153 J[rmap.son2eval_index] = rmap.son_mapping.J[rmap.sons]
154 self.J = J
155 #print(np.min(J), np.max(J))
156 conductivities = ThermicSolver.setconductivity_fromls(conductivity, cut_mesh, ls_tri=ls[self.tris], ls_tri_cut= cut_ls[cut_mesh.new_tris], rmap= rmap)
157 conductivities = np.asarray(conductivities, npfloat)
158 ws = np.asarray(np.eye(2).flatten()*integrationRule0.w[0], npfloat)
159 tmp = np.tensordot(conductivities, ws, axes = 0 )
160 Dval = (J[:, np.newaxis]*tmp).flatten()
161 Dj = np.repeat(np.arange(2*neval).reshape((-1,2)), 2, axis = 0).flatten()
162 Di = np.arange(2*neval).repeat(2)
163 self.D = sp.sparse.csr_array((Dval, (Di, Dj)), (2*neval, 2*neval) )
164 self.Fint = integrationRule0.w[0]*source *self.N.T.dot(J)
165 self.rmap = rmap
166 else :
167 self.tris_xy = self.xy[self.tris]
168 self.tris_map = mapping.MappingT3(self.tris_xy)
169 if space =='P1C':
170 self.space = fs.FEMSpaceP1("Temperature", sm.SMesh(self.xy, self.tris), vecdim = 1, tris_map = self.tris_map)
171 elif space =='P1NC':
172 self.space = fs.FEMSpaceP1NC("Temperature", sm.SMesh(self.xy, self.tris), vecdim = 1, tris_map = self.tris_map)
173 else :
174 print('Error : space ', space, ' not coded in ', ThermicSolver.__name__,'.' )
175 raise NotImplementedError(f"space not coded")
176 self.B = self.space.operator_dof2grad_tri(integrationRule0.uv)
177 self.B0 = self.B
178 nt = len(self.tris)
179 J = self.tris_map.J
180 self.J = J
181 Dval = np.tensordot(J, np.eye(2).flatten()*conductivity*integrationRule0.w[0], axes = 0).flatten()
182 Dj = np.repeat(np.arange(2*nt).reshape((-1,2)), 2, axis = 0).flatten()
183 Di = np.arange(2*nt).repeat(2)
184 self.D = sp.sparse.csr_array((Dval, (Di, Dj)), (2*nt, 2*nt) )
185 self.N = self.space.operator_dof2val_tri(integrationRule0.uv)
186 self.NatVertices = self.space.operator_dof2val_tri(quad.T3_nodes().uv)
189 if callable(source) :
190 quadrature = quad.T3_gauss(5)
191 npt = len(quadrature.w)
192 uv = quadrature.uv
193 weight = quadrature.w
194 NP1 = fs.FEMSpaceP1.N(uv)
195 tris_gp_xy = np.einsum('ijk, ikl-> ijl', NP1[np.newaxis, :,:], self.tris_xy)
196 sourcegp = source(tris_gp_xy)
197 wsourcegp = weight[np.newaxis,:] *J[:, np.newaxis]* sourcegp
198 Nsource = self.space.operator_dof2val_tri(uv)
199 self.Fint = Nsource.T.dot(wsourcegp.flatten())
200 elif np.isscalar(source):
201 self.Fint = integrationRule0.w[0]*source *self.N.T.dot(J)
202 else :
203 raise ValueError("the source term must be a scalar or be evaluatable")
204 def solve(self, solver = "default", compute_rcond = False):
205 if type(self.space) == fs.SumSpace : space = self.space.list_of_space[0]
206 else : space = self.space
207 self.dirichlet2dofs, self.dirichlet2xy = getDofsAndInterpolationPoints(self.modelName, self.physical_dict, self.dirichlet.keys(), space)
208 fixeddofs = np.hstack([ np.array(dofs, dtype = np.int64) for dofs in self.dirichlet2dofs.values() ])
209 fixedvalues = []
210 for xy, dofs, dval in zip(self.dirichlet2xy.values(), self.dirichlet2dofs.values(), self.dirichlet.values()):
211 if np.isscalar(dval) :
212 dval = cast(npfloat, dval)
213 fixedvalues.append(np.ones(len(dofs))*dval)
214 elif callable(dval) :
215 fixedvalues.append(dval(xy))
216 else :
217 raise ValueError("val mus be a scalar or a callable")
218 fixedvalues = np.hstack(fixedvalues)
219 freedof = np.setdiff1d(np.arange(self.B.shape[1]), fixeddofs)
220 Bfree = self.B[:, freedof]
221 U = np.zeros(self.B.shape[1])
222 U[fixeddofs] = fixedvalues
223 gradTdiri = self.B[:, fixeddofs].dot(fixedvalues)
224 Fdiri = - Bfree.T.dot(self.D.dot(gradTdiri))
225 F = self.Fint[freedof] + Fdiri
226 if solver == "default" :
227 Kfree = Bfree.T.dot(self.D).dot(Bfree)
228 luKfree = sp.sparse.linalg.splu(Kfree)
229 U[freedof] = luKfree.solve(F)
230 if compute_rcond :
231 rcond = st.rcond_estimate(Kfree, luKfree)
232 return U, rcond
233 return U
235 if solver == "cholmod":
236 Kfree = Bfree.T.dot(self.D).dot(Bfree)
237 chol = sksparse.cholmod.cholesky(Kfree)
238 U[freedof] = chol(F)
239 if compute_rcond:
240 diag = chol.D()
241 rcond = (np.min(diag)/np.max(diag))**2
242 return U, rcond
243 return U
244 if solver == "AATcholmod":
245 sqrtD = sp.sparse.csr_array( (np.sqrt(self.D.data), self.D.indices, self.D.indptr), self.D.shape)
246 chol = sksparse.cholmod.cholesky_AAt( (sqrtD.dot(Bfree)).T)
247 U[freedof] = chol(F)
248 if compute_rcond :
249 raise NotImplementedError('rcond is not implemeted for AATcholmod')
250 return U
252 print('solver '+ solver + ' is not known !')
253 raise ValueError('solver can only be one of "default", "cholmod", "AATcholmod"' )
255 def postproVal(self, U):
256 return self.tris_xy, self.NatVertices.dot(U).reshape((-1, 3))
257 def postproGrad(self, U):
258 return self.tris_xy, self.B.dot(U).reshape((-1, 2))
259 def H1_error(self, U, gradexactfun, quadrature = quad.T3_gauss(0)):
260 npt = len(quadrature.w)
261 uv = quadrature.uv
262 weight = quadrature.w
263 B = self.space.operator_dof2grad_tri(uv, self.rmap)
264 gradNum = B.dot(U).reshape((-1, npt, 2))
265 NP1 = fs.FEMSpaceP1.N(uv)
266 tris_gp_xy = np.einsum('ijk, ikl-> ijl', NP1[np.newaxis, :,:], self.tris_xy)
267 gradExact = gradexactfun(tris_gp_xy)
268 error2_gp = np.linalg.norm((gradNum - gradExact), axis = 2)**2
269 error2 = np.sum(self.J*np.einsum( 'ij, ij-> i', weight[np.newaxis, :] , error2_gp))
270 return np.sqrt(error2)
271 def maxValError(self, U, exactfun, samplingpoints = quad.T3_gauss(0).uv) :
272 uv = samplingpoints
273 npt = len(uv)
274 N = self.space.operator_dof2val_tri(uv, self.rmap)
275 valNum = N.dot(U).reshape((-1, npt))
276 NP1 = fs.FEMSpaceP1.N(uv)
277 tris_gp_xy = np.einsum('ijk, ikl-> ijl', NP1[np.newaxis, :,:], self.tris_xy)
278 valExact = exactfun(tris_gp_xy)
279 error_gp = np.abs(valNum - valExact)
280 return np.max(error_gp)
281 def maxGradError(self, U, gradexactfun, samplingpoints = quad.T3_gauss(0).uv) :
282 uv = samplingpoints
283 npt = len(uv)
284 B = self.space.operator_dof2grad_tri(uv, self.rmap)
285 gradNum = B.dot(U).reshape((-1, npt, 2))
286 NP1 = fs.FEMSpaceP1.N(uv)
287 tris_gp_xy = np.einsum('ijk, ikl-> ijl', NP1[np.newaxis, :,:], self.tris_xy)
288 gradExact = gradexactfun(tris_gp_xy)
289 error_gp = np.linalg.norm((gradNum - gradExact), axis = 2)
290 return np.max(error_gp)
291 def L2_error(self, U, exactfun, quadrature = quad.T3_gauss(0)):
292 npt = len(quadrature.w)
293 uv = quadrature.uv
294 weight = quadrature.w
295 N = self.space.operator_dof2val_tri(uv, self.rmap)
296 ValNum = N.dot(U).reshape((-1, npt))
297 NP1 = fs.FEMSpaceP1.N(uv)
298 tris_gp_xy = np.einsum('ijk, ikl-> ijl', NP1[np.newaxis, :,:], self.tris_xy)
299 ValExact = exactfun(tris_gp_xy)
300 error2_gp = (ValNum - ValExact)**2
301 error2 = np.sum(self.J*np.einsum( 'ij, ij-> i', weight[np.newaxis, :] , error2_gp))
302 return np.sqrt(error2)
304if __name__== '__main__':
305 import gmsh
306 from pyxfem.core.exec_tools import Profile
307 import pyxfem.core.gmsh_mesh_generator as mg
308 import pyxfem.postpro.gmsh_post_pro_interface as gp
310 elementSize = 0.2
311 solver = "default" #"cholmod" # "default", "AATcholmod"
312 popUpGmsh = True
313 Profile.doprofile = False
315 dofem = ['P1C', 'P1NC']
316 dofem2conductivity = True
317 doxfem2conductivity = ['xFEMSpace', 'xRidgeSpaceP1']
318 doxfem_bench = True
320# dofem = ['P1C']
321# dofem2conductivity = False
322# doxfem2conductivity = []
323# doxfem_bench = False
325 gmsh.initialize()
326 gmsh.option.setNumber('PostProcessing.AnimationCycle', 1)
327 gmsh.option.setNumber('Mesh.SurfaceEdges', 1 if elementSize > 0.1 else 0)
328 modelName, physical_dict = mg.l_shape(elementSize)
329 nv = mi.nb_vertices(modelName)
330 nf = mi.nb_triangles(modelName)
331 print('nb elements =', nf)
332 print('nb vertices =', nv)
333 physname_to_dirichlet_value = {'dOmega_left':0.33, 'dOmega_right' :0.1}
334 conductivity = 1.
335 for spaceName in dofem :
336 print('Solving a fem problem using space ', spaceName )
337 pbfem = ThermicSolver(modelName, physical_dict, conductivity, physname_to_dirichlet_value, space = spaceName, enrich = False )
338 U = pbfem.solve(solver = solver)
339 trisxy, Ttris = pbfem.postproVal(U)
340 trisxy, GTtris = pbfem.postproGrad(U)
341 gp.listPlotFieldTri(trisxy, Ttris, P0 = False, viewname ='T_'+ spaceName, NbIso = 20)
342 gp.listPlotFieldTri(trisxy, GTtris, P0 = True, viewname ='GradT_'+ spaceName, NbIso = 20, VectorType = "Arrow")
343 space = pbfem.space
344 if type(space) == fs.FEMSpaceP1NC:
345 N_edge = space.operator_edgetrace(quad.Edge_nodes())
346 valEdge = N_edge.dot(U).reshape(-1,2)
347 gp.listPlotFieldLine(space.m.e2v_coord(), valEdge[:,0], P0 = False, viewname ='T_EdgeL', LineType = "3D cylinder", LineWidth = 4)
348 gp.listPlotFieldLine(space.m.e2v_coord(), valEdge[:,1], P0 = False, viewname ='T_EdgeR', LineType = "3D cylinder", LineWidth = 4)
349 if dofem2conductivity :
350 print('Solving a fem problem with 2 conductivities defined by a levelset')
351 lsfun = lscut.SimpleLevelSetFunction.disk([-0.7,0.5], 1.)
352 pbfem = ThermicSolver(modelName, physical_dict, [20*conductivity, conductivity], physname_to_dirichlet_value, lsfunforconductivity = lsfun, enrich = False )
353 U = pbfem.solve(solver = solver)
354 trisxy, Ttris = pbfem.postproVal(U)
355 trisxy, GTtris = pbfem.postproGrad(U)
356 gp.listPlotFieldTri(trisxy, Ttris, P0 = False, viewname ='T2cond_fem', NbIso = 20)
357 gp.listPlotFieldTri(trisxy, GTtris, P0 = True, viewname ='GradT2cond_fem', NbIso = 20, VectorType = "Displacement")
358 for enrichmentName in doxfem2conductivity :
359 print('Solving a xfem problem with 2 conductivities defined by a levelset. enrichment is ', enrichmentName)
360 lsfun = lscut.SimpleLevelSetFunction.disk([-0.7,0.5], 1.)
361 pbxfem = ThermicSolver(modelName, physical_dict, [20*conductivity, conductivity], physname_to_dirichlet_value, lsfunforconductivity = lsfun, enrich = enrichmentName)
362 print(" nb enriched dofs = ", pbxfem.xspace.scalar_size())
363 U = pbxfem.solve(solver = solver)
364 trisxy, Ttris = pbxfem.postproVal(U)
365 trisxy, GTtris = pbxfem.postproGrad(U)
366 gp.listPlotFieldTri(trisxy, Ttris, P0 = False, viewname ='T2cond_'+enrichmentName, NbIso = 20)
367 gp.listPlotFieldTri(trisxy, GTtris, P0 = True, viewname ='GradT2cond_'+enrichmentName, NbIso = 20, VectorType = "Displacement")
368 if doxfem_bench :
369 print('Solving a xfem problem with 2 conductivities defined by a more complicated levelset.')
370 x = np.linspace(-0.8, 0.8, 4)
371 y = np.linspace(-0.8, 0.8, 4)
372 lsfun = lscut.SimpleLevelSetFunction.union([lscut.SimpleLevelSetFunction.disk([xi,yi], 0.10) for xi in x for yi in y] )
373 pbxfem = ThermicSolver(modelName, physical_dict, [conductivity*0.001, conductivity], physname_to_dirichlet_value, lsfunforconductivity = lsfun, enrich = 'xFEMSpace')
374 print(" nb enriched dofs = ", pbxfem.xspace.scalar_size())
375 U = pbxfem.solve(solver = solver)
376 trisxy, Ttris = pbxfem.postproVal(U)
377 trisxy, GTtris = pbxfem.postproGrad(U)
378 gp.listPlotFieldTri(trisxy, Ttris, P0 = False, viewname ='BenchXfem_T', NbIso = 20)
379 gp.listPlotFieldTri(trisxy, GTtris, P0 = True, viewname ='BenchXfem_GradT', NbIso = 20, VectorType = "Displacement")
380 if popUpGmsh : gmsh.fltk.run()
381 gmsh.finalize()
382 Profile.print_stats()