Coverage for pyxfem/core/quadrature.py: 28%
117 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 Mon Jan 8 13:43:53 2024
6@author: nchevaug
7"""
8import numpy as np
10class Edge_gauss():
11 def __init__(self, order=0):
12 npt = np.int64(np.ceil((order+1)/2))
13 s, w = np.polynomial.legendre.leggauss(npt)
14 self.s = (1.+s)/2.
15 self.w = w/2.
16class Edge_nodes():
17 def __init__(self, order=0):
18 self.s = np.array([0.,1.])
19 self.w = np.array([0.5,0.5])
20class T3_edges():
21 def __init__(self, edgeRule):
22 s = edgeRule.s
23 npt = len(s)
24 u = np.hstack((s, 1-s, np.zeros(npt) ))
25 v = np.hstack((np.zeros(npt), s, 1-s))
26 self.uv = np.column_stack((u,v))
27 self.w = np.tile(edgeRule.w/3.,3)
28class T3_gauss():
29 _maxorder = 8
30 def __init__(self, order=0, negweightOK = False):
32 if order in [0,1] :
33 self.w = np.array([0.5])
34 self.uv = np.array([[1./3., 1./3.]])
35 return
36 if order == 2 :
37 self.w = np.repeat(1./6.,3)
38 self.uv = np.array([[1./6., 1./6.], [2./3., 1./6.], [1./6, 2./3.] ])
39 return
40 if order == 3 and negweightOK :
41 self.w = np.array([-9./32., 25./96., 25./96, 25/96. ])
42 self.uv = np.array([[1./3., 1./3.], [0.6, 0.2], [0.2,0.6], [0.2, 0.2] ])
43 return
44 if order in [3, 4]:
45 self.w = np.repeat([ 0.109951743655322/2., 0.223381589678011/2.], 3)
46 self.uv = np.array([
47 [0.816847572980459, 0.091576213509771], [0.091576213509771, 0.816847572980459],
48 [0.091576213509771, 0.091576213509771], [0.108103018168070, 0.445948490915965],
49 [0.445948490915965, 0.108103018168070], [0.445948490915965, 0.445948490915965]])
50 return
51 if order == 5 :
52 self.w = 0.5*np.array([0.225, 0.125939180544827, 0.125939180544827, 0.125939180544827, 0.132394152788506, 0.132394152788506, 0.132394152788506])
53 self.uv = np.array([
54 [0.333333333333333, 0.333333333333333], [0.797426985353087, 0.101286507323456],
55 [0.101286507323456, 0.797426985353087], [0.101286507323456, 0.101286507323456],
56 [0.470142064105115, 0.059715871789770], [0.059715871789770, 0.470142064105115],
57 [0.470142064105115, 0.470142064105115]])
58 return
59 if order == 6 :
60 self.w = 0.5*np.array([
61 0.050844906370207, 0.050844906370207, 0.050844906370207,
62 0.116786275726379, 0.116786275726379, 0.116786275726379,
63 0.082851075618374, 0.082851075618374, 0.082851075618374,
64 0.082851075618374, 0.082851075618374, 0.082851075618374])
65 self.uv = np.array([
66 [0.873821971016996, 0.063089014491502], [0.063089014491502, 0.873821971016996],
67 [0.063089014491502, 0.063089014491502], [0.501426509658179, 0.249286745170910],
68 [0.249286745170910, 0.501426509658179], [0.249286745170910, 0.249286745170910],
69 [0.636502499121399, 0.310352451033785], [0.310352451033785, 0.636502499121399],
70 [0.636502499121399, 0.053145049844816], [0.310352451033785, 0.053145049844816],
71 [0.053145049844816, 0.310352451033785], [0.053145049844816, 0.636502499121399]])
72 return
73 if order == 7 and negweightOK :
74 self.w = np.array([
75 -0.149570044467682, 0.175615257433208, 0.175615257433208,
76 0.053347235608838, 0.053347235608838, 0.053347235608838,
77 0.077113760890257, 0.077113760890257, 0.077113760890257,
78 0.077113760890257, 0.077113760890257, 0.077113760890257])
79 self.uv = np.array([
80 [0.333333333333333, 0.333333333333333],
81 [0.479308067841920, 0.260345966079040],
82 [0.260345966079040, 0.479308067841920],
83 [0.869739794195568, 0.065130102902216],
84 [0.065130102902216, 0.869739794195568],
85 [0.065130102902216, 0.065130102902216],
86 [0.048690315425316, 0.312865496004874],
87 [0.312865496004874, 0.048690315425316],
88 [0.638444188569810, 0.048690315425316],
89 [0.048690315425316, 0.638444188569810],
90 [0.312865496004874, 0.638444188569810],
91 [0.638444188569810, 0.312865496004874],
92 ])
93 return
94 if order in [7,8] :
95 self.w = np.array([
96 0.144315607677787, 0.095091634267285, 0.095091634267285,
97 0.095091634267285, 0.103217370534718, 0.103217370534718,
98 0.103217370534718, 0.032458497623198, 0.032458497623198,
99 0.032458497623198, 0.027230314174435, 0.027230314174435,
100 0.027230314174435, 0.027230314174435, 0.027230314174435,
101 0.027230314174435])
102 self.uv = np.array([
103 [0.333333333333333, 0.333333333333333],
104 [0.081414823414554, 0.459292588292723],
105 [0.459292588292723, 0.081414823414554],
106 [0.459292588292723, 0.459292588292723],
107 [0.658861384496480, 0.170569307751760],
108 [0.170569307751760, 0.658861384496480],
109 [0.170569307751760, 0.170569307751760],
110 [0.898905543365938, 0.050547228317031],
111 [0.050547228317031, 0.898905543365938],
112 [0.050547228317031, 0.050547228317031],
113 [0.008394777409958, 0.728492392955404],
114 [0.728492392955404, 0.008394777409958],
115 [0.263112829634638, 0.008394777409958],
116 [0.008394777409958, 0.263112829634638],
117 [0.263112829634638, 0.728492392955404],
118 [0.728492392955404, 0.263112829634638]])
119 return
120 raise Exception('T3_gausspoints are only defined up to order {:d}'.format(T3_gauss._maxorder))
121class T3_nodes():
122 def __init__(self):
123 self.w = np.repeat(0.5/3.,3)
124 self.uv = np.array([[0.,0.], [1.,0.], [0.,1.]])
125class T3_midedges():
126 def __init__(self):
127 self.w = np.repeat(1./6.,3)
128 self.uv = np.array([[0.5, 0.], [0.5,0.5], [0.,0.5]])
130if __name__== '__main__' :
131 import gmsh
132 import pyxfem.gmsh_post_pro_interface as gp
134 T0 = np.array([[0.,0.], [1.,0.], [0.,1.]])
135 T1 = np.array([[1.,0.], [1.,1.], [0.,1.]])
136 J = np.array([1.,1.])
137 tris = np.array([T0, T1])
138 fun = lambda xy : np.ones(xy.shape[:-1])
139 fun = lambda xy : xy[...,0]**5
140 for i in np.arange(6) :
141 quad = T3_gauss(i)
142 u, v = quad.uv[:,0], quad.uv[:,1]
143 N = np.array([(1.-u-v), u, v])
144 egxy = np.einsum ('esg, esd -> egd', N[np.newaxis,:,:], tris)
145 egval = fun(egxy)
146 integral = np.einsum('eg, eg', J[:, np.newaxis].dot(quad.w[np.newaxis,:]), egval )
147 print('integration order :', i, 'integrale', integral)
151 popUpGmsh = True
152 gmsh.initialize()
153 gmsh.fltk.initialize()
155 xy = np.array([[0.,0.],[1.,0.],[0.,1.]])
157 edge = np.array([[0,1]])
158 xyedge = xy[edge]
159 for order in range(9) :
160 rule = Edge_gauss(order)
161 npt = len(rule.s)
162 s = rule.s.reshape((npt,1))
163 w = rule.w.reshape((npt,1))
164 view = gp.listPlotFieldPoint(s, w, viewname= 'Order = '+ str(order) , PointType='Scaled sphere', PointSize= 10)
165 gp.listPlotFieldLine(xyedge, np.array([0.001]), P0 = True, LineType='3D cylinder', LineWidth=1, gv = view)
167 tri = np.array([[0,1,2]])
168 xyTri = xy[tri]
169 edges = np.array([[0,1],[1,2],[2,0]])
170 xyEdges = xy[edges]
172 for order in range(T3_gauss._maxorder+1) :
173 print(order)
174 rule = T3_gauss(order)
175 npt = len(rule.uv)
176 uv = rule.uv.reshape((npt,2))
177 w = rule.w.reshape((npt,1))
178 view = gp.listPlotFieldPoint(uv, w, viewname= 'Order = '+ str(order) , PointType='Scaled sphere', PointSize= 10)
179 gp.listPlotFieldLine(xyEdges, np.array([0.001]*3), P0 = True, LineType='3D cylinder', LineWidth=1, gv = view)
181 for order in range(T3_gauss._maxorder+1) :
182 rule = T3_edges(Edge_gauss(order))
183 npt = len(rule.uv)
184 uv = rule.uv.reshape((npt,2))
185 w = rule.w.reshape((npt,1))
186 view = gp.listPlotFieldPoint(uv, w, viewname= 'Order = '+ str(order) , PointType='3D sphere', PointSize= 10)
187 gp.listPlotFieldLine(xyEdges, np.array([0.001]*3), P0 = True, LineType='3D cylinder', LineWidth=1, gv = view)
189 gmsh.option.setNumber('PostProcessing.AnimationCycle', 1)
190 #gmsh.open('post/quadrature.opt')
191 #gp.gmshExportsViews('post/quadrature', 'mpg')
192 #gp.gmshExportsViews('post/quadrature', 'pos')
193 if popUpGmsh:
194 gmsh.fltk.run()
195 #gp.gmshExportsViews('post/quadrature', 'jpg')
196 #gmsh.write('post/quadrature.opt')
198 gmsh.fltk.finalize()
199 gmsh.finalize()