Coverage for pyxfem/core/mapping.py: 84%
102 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 Dec 18 14:05:29 2023
6@author: nchevaug
7"""
8import numpy as np
10from pyxfem.core.config import IndexType
12def tri2g(t3, dim = None):
13 ''' transform the node coordinate of an array of triangles [..., 3,2|3]
14 to an array g[...,2,2|3], where g[...,i, :] is the vector linking vertice 0 to vertice i
15 for each each entry in the array ...'''
16 assert t3.shape[-2] == 3
17 if dim is None:
18 dim =t3.shape[-1]
19 return (t3[...,1:, :dim] - t3[...,0,np.newaxis,:dim]).reshape(t3.shape[:-2] + (2,dim))
21class MappingT3():
22 ''' construct mapping from an array of T3 nodes coordinates '''
23 def __init__(self, t3xy, hmin=0., alpha=0.1, p=1):
24 # pylint: disable=C0103
25 self.hmin = hmin
26 self.alpha = alpha
27 self.p = p
28 self.Jmin = alpha*((hmin)**p)**2
29 self.xy = t3xy
30 #self.F = np.swapaxes(self.xy[..., 1:, :]-self.xy[..., [0], :], -1, -2)
31 self.F = tri2g(t3xy, 2).swapaxes(1,2)
32 F = self.F
33 self.J = F[:, 0, 0] * F[:, 1, 1] - F[:, 1, 0] * F[:, 0, 1]
34 limit_reached = self.J < self.Jmin
35 self.Jlimited = np.where(limit_reached, self.Jmin, self.J)
36 self.cofac_F = np.swapaxes(
37 np.array([[F[:, 1, 1], -F[:, 1, 0]], [-F[:, 0, 1], F[:, 0, 0]]]), 0, 2)
38 self.inv_F = self.cofac_F / self.J[:, np.newaxis, np.newaxis]
39 self.inv_F_limited = self.cofac_F/self.Jlimited[:, np.newaxis, np.newaxis]
41 def uv2xy(self, uv, elem_ids=None):
42 """ map uv local coordinates to xy for all triangles """
43 if elem_ids is None:
44 return self.xy[..., [0], :] + np.einsum('...jl,...kl->...kj', self.F, uv)
45 else:
46 return self.xy[elem_ids][..., [0], :] \
47 + np.einsum('...jl,...kl->...kj', self.F[elem_ids], uv)
49 def xy2uv(self, xy, elem_ids=None):
50 """ return uv in each triangle for each entry of xy """
51 if elem_ids is not None:
52 invFs = self.inv_F[elem_ids] # pylint: disable=C0103
53 xy0s = xy - self.xy[elem_ids][..., [0], :]
54 else:
55 invFs = self.inv_F # pylint: disable=C0103
56 xy0s = xy - self.xy[..., [0], :]
57 return np.einsum('...jl,...kl->...kj', invFs, xy0s)
60class RemapHelper():
61 ''' a class to help mapping data when we have a partition of some 'parent element'
62 into a set of sons element'''
64 def __init__(self, index_parent2sons, sons, parent_mapping, son_mapping):
65 self.index_parent2sons = index_parent2sons
66 self.sons = sons
67 self.parent_mapping = parent_mapping
68 self.son_mapping = son_mapping
70 self.nbappro = len(index_parent2sons)-1
71 self.indextype = index_parent2sons.dtype
72 self.appro2nbsons = index_parent2sons[1:] - index_parent2sons[:-1] # pylint: disable=C0103
73 # elements for which eval element is the appro element (appro as no sons)
74 self.approEQeval = np.where(self.appro2nbsons == 0)[0] # pylint: disable=C0103
75 # elements for which eval element is the appro element (appro as no sons)
76 self.approNEQeval = np.where(self.appro2nbsons > 0)[0] # pylint: disable=C0103
77 self.appro2nbevals = np.where(
78 self.appro2nbsons == 0, 1, self.appro2nbsons)
79 appro2evals_index = np.zeros(len(self.appro2nbevals)+1, dtype= IndexType)
80 np.cumsum(self.appro2nbevals, dtype= IndexType, out=appro2evals_index[1:])
81 self.nbevals = appro2evals_index[-1]
82 self.noson2eval_index = appro2evals_index[self.approEQeval]
83 self.son2eval_index = np.setdiff1d(np.arange(self.nbevals), self.noson2eval_index)
84 self.eval2appro = np.repeat( np.arange(self.nbappro), self.appro2nbevals)
85 self.son2parent_elem_ids = np.repeat(self.approNEQeval, self.appro2nbsons[self.approNEQeval])
87 def remap_on_parent(self, evaltrisuv):
88 """ compute local coordinate on parent element from local coordinate on son elements """
89 sons_point_xy = self.son_mapping.uv2xy(evaltrisuv)[self.sons]
90 sons_point_on_parents_uv = self.parent_mapping.xy2uv(
91 sons_point_xy, elem_ids=self.son2parent_elem_ids)
92 return sons_point_on_parents_uv
94 def get_eval_tris(self):
95 """ retun an np.array containing the coordinates of the nodes of each sons elements
96 if an element has no sons, it's the coordinates of this element that are returned
97 The elements are sorted according to there parent element index
98 """
99 xytris = np.zeros((self.nbevals,) + (3, 2))
100 xytris[self.noson2eval_index] = self.parent_mapping.xy[self.approEQeval]
101 xytris[self.son2eval_index] = self.son_mapping.xy[self.sons]
102 return xytris
105class RemapHelperFiltered():
106 """ a class to help mapping data when we have a partition of some 'parent element'
107 into a set of sons element.
108 This version include a filter on parent element, only on parent
109 """
111 def __init__(self, index_parent2sons, sons, parent_mapping, son_mapping, only_on_parents):
112 self.indextype = index_parent2sons.dtype
113 self.nbappro = len(only_on_parents)
114 self.only_on_parents = only_on_parents
115 self.sons = np.concatenate(np.split(sons, np.column_stack(
116 (index_parent2sons[only_on_parents],
117 index_parent2sons[only_on_parents+1])).flatten())[1::2])
118 self.appro2nbsons = index_parent2sons[only_on_parents +1] \
119 - index_parent2sons[only_on_parents]
120 self.parent_mapping = parent_mapping
121 self.son_mapping = son_mapping
122 self.approEQeval_all = np.where(self.appro2nbsons == 0)[0] # pylint: disable=C0103
123 self.approNEQeval_all = np.where(self.appro2nbsons > 0)[0] # pylint: disable=C0103
125 # elements for which eval element is the appro element (appro as no sons)
126 self.approEQeval = only_on_parents[self.approEQeval_all] # pylint: disable=C0103
127 # elements for which eval element is the appro element (appro as no sons)
128 self.approNEQeval = only_on_parents[self.approNEQeval_all] # pylint: disable=C0103
129 self.appro2nbevals = np.where(self.appro2nbsons == 0, 1, self.appro2nbsons)
130 appro2evals_index = np.zeros(len(self.appro2nbevals)+1, dtype=IndexType)
131 np.cumsum(self.appro2nbevals, dtype=IndexType, out=appro2evals_index[1:])
132 self.nbevals = appro2evals_index[-1]
133 self.noson2eval_index = appro2evals_index[self.approEQeval_all]
134 self.son2eval_index = np.setdiff1d(np.arange(self.nbevals), self.noson2eval_index)
135 self.eval2appro = np.repeat(np.arange(self.nbappro), self.appro2nbevals)
136 self.son2parent_elem_ids = self.only_on_parents[np.repeat(
137 self.approNEQeval_all, self.appro2nbsons[self.approNEQeval_all])]
139 def remap_on_parent(self, evaltrisuv):
140 """ compute local coordinate on parent element from local coordinate on son elements """
141 sons_point_xy = self.son_mapping.uv2xy(evaltrisuv, self.sons)
142 sons_point_on_parents_uv = self.parent_mapping.xy2uv(
143 sons_point_xy, elem_ids=self.son2parent_elem_ids)
144 return sons_point_on_parents_uv
146 def get_eval_tris(self):
147 """ retun an np.array containing the coordinates of the nodes of each sons elements
148 if an element has no sons, it's the coordinates of this element that are returned
149 The elements are sorted according to there parent element index
150 """
151 xytris = np.zeros((self.nbevals,) + (3, 2))
152 xytris[self.noson2eval_index] = self.parent_mapping.xy[self.approEQeval_all]
153 xytris[self.son2eval_index] = self.son_mapping.xy[self.sons]
154 return xytris
158class Mapping3DTri():
159 ''' mapping of 3d triangles with 3d coordinates.'''
160 def __init__(self, t3):
161 assert(t3.shape[-2:] in [(3,2), (3,3)])
162 # pylint: disable=C0103
163 self.g = tri2g(t3)
164 self.m = self.g@self.g.swapaxes(-1,-2)
165 self.Q =np.zeros( t3.shape[:-2] + (3,3))
166 self.Q[...,:2,:t3.shape[-1]] = self.g
167 self.Q[..., 2,:t3.shape[-1]] = np.cross(self.Q[...,0,:], self.Q[...,1,:])
168 self.Q[..., 1,:t3.shape[-1]] = np.cross(self.Q[...,2,:], self.Q[...,0,:])
169 self.Q /= np.linalg.norm(self.Q, axis = -1, keepdims=True)
170 self.F2D = self.Q[..., :-1, :]@self.g.swapaxes(-1,-2)
172# if __name__ == '__main__':
173# #pass
174# # tri_xy = np.array([[0.1,0.3],[1.,0.],[0.5,0.5]])
175# # print('triangle :')
176# # print(tri_xy)
177# # t3m = t3map_old(tri_xy)
178# # print('_____ direct, using map ')
179# # print('area =', np.cross(tri_xy[1] - tri_xy[0], tri_xy[2] - tri_xy[0])/2., t3m.detJ()/2. )
180# # print('F', (tri_xy[1:,:] - tri_xy[0,:]).T, t3m.F())
181# # uv_list = [ np.array([[1./3., 1./3.]]),
182# # np.array([[1./3., 1./3.], [0.,0.], [1.,0.], [0.,1.]])]
183# # for i, uv in enumerate(uv_list) :
184# # print('test {}'.format(i))
185# # print(' uv')
186# # print(uv)
187# # xy = t3m.uv2xy(uv)
188# # uv = t3m.xy2uv(xy)
189# # print(' xy')
190# # print( xy)
191# # print(' uv')
192# # print( uv)
193# #
194# triangles = np.array([[[0.,0., 2.], [1.,0.,2.], [0,1,2]] , [[0.,0.,0.],[1.,0.,0.], [0,1,0]]])
195# #t3 = np.array([[0.1,0.2,0.12], [1.,0.,0.], [0,1,0]])
196# #map = mapping3dtri(tri_xyz)
197# g = tri2g(triangles)
198# mapp = Mapping3DTri(triangles)
199# print(mapp.g)
200# print(mapp.Q)
201# print(mapp.F2D)