pyviewfactor.pyviewfactor

  1""" PVF main functions"""
  2
  3
  4# Imports
  5from numba import njit
  6import pyvista as pv
  7import pandas as pd
  8import scipy.integrate
  9import numpy as np
 10
 11
 12# Functions
 13def fc_unstruc2poly(mesh_unstruc):
 14    """Convenience conversion function from UnstructuredGrid to PolyData
 15
 16    Parameters
 17    ----------
 18    * **mesh_unstruc** : *pyvista.UnstructuredGrid*
 19
 20        > Unstructured Pyvista Grid.
 21
 22    Returns
 23    -------
 24    * *pyvista.PolyData*
 25
 26        > The same mesh converted to a surface pyvista.PolyData.
 27
 28    Examples
 29    --------
 30    >>> import pyviewfactor as pvf
 31    >>> import pyvista as pv
 32    >>> sphere = pv.Sphere(radius=0.5, center=(0, 0, 0))
 33    >>> subset = sphere.extract_cells(10)
 34    >>> subsetPoly = fc_unstruc2poly(subset)
 35    >>> subsetPoly
 36    PolyData (0x1fdd9786040)
 37      N Cells:    1
 38      N Points:    3
 39      X Bounds:    -5.551e-17, 3.617e-02
 40      Y Bounds:    0.000e+00, 4.682e-02
 41      Z Bounds:    -5.000e-01, -4.971e-01
 42      N Arrays:    0
 43
 44    """
 45
 46    # Get the points and cells
 47    points = mesh_unstruc.points
 48    faces = mesh_unstruc.cells
 49    # Return the same geometry as a pv.PolyData mesh
 50    return pv.PolyData(points, faces)
 51
 52
 53def get_visibility(c1, c2, strict=False,print_warning=True,
 54                         rounding_decimal = 1):
 55    """Facets visibility:
 56
 57    A test to check if two facets can "see" each other, taking the normals
 58    into consideration (no obstruction tests, only normals orientations).
 59
 60    Parameters
 61    ----------
 62    * **c1** : *pyvista.PolyData*
 63
 64        > PolyData facet (pyvista format).
 65
 66    * **c2** : *pyvista.PolyData*
 67
 68        > PolyData facet (pyvista format).
 69
 70    * **strict** : *Bool*
 71
 72        > If *True*, checks all the points are able to see each other
 73          (considering the face normal) and then continue. If some points are
 74          "begind the other faces, it will return *False*,
 75        > Else, compute the centroids visibility, and might print a warning if
 76          some points are "behind".
 77
 78    * **print_warning** : *Bool*
 79
 80        > If *True*, warning messages will be printed in addition to be returned
 81
 82    * **rounding_decimal** : *Int*
 83
 84        > Number of decimals at which rounding shall be done for the points
 85          coordinates. This is to avoid numeric issues when the points
 86          coordinates are given as for instance 1.234e-13 instead of 0.
 87
 88    Returns
 89    -------
 90    * **bool**
 91
 92        > True when the facets "see" each other, False else.
 93
 94    * **str**
 95
 96        > Warning message if any (empty string if no warning).
 97
 98    Examples
 99    --------
100    >>> import pyvista as pv
101    >>> import pyviewfactor as pvf
102    >>> tri1 = pv.Triangle([[0.0, 1.0, 1.0],[1.0, 1.0, 1.0],[1.0, 0.0, 1.0]])
103    >>> tri2 = pv.Triangle([[1.0, 0.0, 0.0],[1.0, 1.0, 0.0],[0.0, 1.0, 0.0]])
104    >>> pvf.get_visibility(tri1, tri2,strict=False, print_warning=True)
105    True
106    """
107
108    # Checking if every point of c1 in "in front" of c2 with the sign of the
109    # triple product:
110    # for pt in c1.points:
111    #     v = pt - c2.points[0]
112    #     edge2a = c2.points[0] - c2.points[1]
113    #     edge2b = c2.points[0] - c2.points[2]
114    #     # "triple" product computation --> "v.(edge2a x edge2b)"
115    #     det = v.dot(np.cross(edge2a, edge2b))
116    #     if det < 0:
117    #         if not strict:
118    #             print(
119    #                 "... ! ... PVF Warning : at least one point of cell is \"behind\" cell2\n")
120    #         else:
121    #             return False
122
123    # Test on the normals orientations for the visibility:
124    # Get cell centers
125    center1 = c1.cell_centers().points
126    center2 = c2.cell_centers().points
127    # Get the vector between the cell centers
128    v21 = center1-center2
129    # Get cell normals
130    n1 = c1.cell_normals
131    n2 = c2.cell_normals
132    # Compute the 2 scalar product
133    pos_dot_prod = np.einsum('ij,ij->i', v21, n2)
134    neg_dot_prod = np.einsum('ij,ij->i', v21, n1)
135    # Return result of visibility test
136    if not (pos_dot_prod > 0 and neg_dot_prod < 0):
137        return (False, "")
138
139    ## Check if the cells are entirely visible to each other
140    ## Checking if every point of c1 in "in front" of c2 with the sign of the
141    ## triple product, and reciprocally
142    for cel_i in [c1, c2]:
143        if cel_i == c1:
144            cel_j = c2
145        else:
146            cel_j = c1
147        edge2a = np.round(cel_j.points[0] - cel_j.points[1], rounding_decimal)
148        edge2b = np.round(cel_j.points[0] - cel_j.points[2], rounding_decimal)
149        edges_cross = np.cross(edge2a, edge2b)
150        det_is_pos = False
151        det_is_neg = False
152        for pt in cel_i.points:
153            v = np.round(pt - cel_j.points[2], rounding_decimal)
154            # "triple" product computation --> "v.(edge2a x edge2b)"
155            det = v.dot(edges_cross)
156            if det < 0:
157                det_is_neg = True
158            elif det > 0:
159                det_is_pos = True
160            if det_is_neg and det_is_pos: # det < 0:
161                if strict:
162                    warning_str = (
163                        "[PVF-Warning-1] strict being True, cells are considered "
164                        "not visible although they partially are"
165                        )
166                    if print_warning:
167                        print(warning_str)
168                    return (False, warning_str)
169                else:
170                    warning_str = (
171                        "[PVF-Warning-2] strict being False, cells are considered "
172                        "entirely visible although they are partially hidden")
173                    if print_warning:
174                        print(warning_str)
175                    return (True, warning_str)
176    return (True, "")
177
178
179def get_visibility_raytrace(face1, face2, obstacle):
180    """Raytrace between face1 and face2
181
182    A test to check if there is an obstruction between two facets.
183
184    Parameters
185    ----------
186    * **face1** : *pyvista.PolyData*
187
188        > face1 to be checked for obstruction.
189
190    * **face2** : *pyvista.PolyData*
191
192        > face2 to be checked for obstruction.
193
194    * **obstacle** : *pyvista.PolyData*
195
196        > The mesh inbetween, composing the potential obstruction.
197
198    Returns
199    -------
200    * *bool*
201
202        > True if no obstruction, False else.
203
204    Examples
205    --------
206    >>> import pyvista as pv
207    >>> import pyviewfactor as pvf
208    >>> tri1 = pv.Triangle([[0.0, 1.0, 1.0],[1.0, 1.0, 1.0],[1.0, 0.0, 1.0]])
209    >>> tri2 = pv.Triangle([[1.0, 0.0, 0.0],[1.0, 1.0, 0.0],[0.0, 1.0, 0.0]])
210    >>> obstacle = pv.Circle(radius=3.0)
211    >>> obstacle.translate([0, 0, 0.5], inplace = True)
212    >>> pvf.get_visibility_raytrace(tri2, tri1, obstacle)
213    False
214
215    """
216    # Define line segment from one cell center to the other
217    start = face1.cell_centers().points[0]
218    stop = face2.cell_centers().points[0]
219    # Perform ray trace along the line segment
220    # points : location of the intersection
221    # ind : indices if the intersection cells
222    points, ind = obstacle.ray_trace(start, stop, first_point=False)
223    # If considering a single cell
224    # /!\ Ce n'est pas le bon test, il faut checker si les points de départ
225    # et d'arrivé appartiennet à "obstacle" !
226    # tester si le obstacle est un maillage ouvert ?
227    # tester
228    if obstacle.n_cells == 1:
229        # The cells are not obstructed if the ray trace from cell 1 to cell2
230        # does *not* intersect the "obstacle" mesh, >> ind is empty
231        return True if ind.size == 0 else False
232    # Il manque un cas de figure
233
234    # Else, if face1 and face2 are contained in the obstacle mesh
235    else:
236        # The cells are obstructed if the ray trace from cell 1 to cell 2 hits
237        # the "obstacle" more than 3 times (on cell 1, cell 2 and at least once
238        # somewhere in between)
239        # If the cells are not obstructed, len(ind) should be == 2
240        return False if len(ind) > 2 else True
241
242
243def get_visibility_MT(face1, face2, obstacle):
244    """Triangle / Ray intersection based on the Möeller-Trumbore algorithm
245
246    A test to check if there is an obstruction between two facets, basaed on
247    Möller-Trumbore algorithm
248
249    **To add**: a "strict" argument, to check not the centroids, but all points
250    of face1 against all point of face2 given an obstacle.
251
252    Parameters
253    ----------
254    * **face1** : *pyvista.PolyData*
255
256        > face1 to be checked for obstruction. *A single PolyData face*
257
258    * **face2** : *pyvista.PolyData*. *A single PolyData face*
259
260        > face2 to be checked for obstruction.
261
262    * **obstacle** : *pyvista.PolyData*. *A single PolyData face, or a entire
263        mesh *
264
265        > The mesh inbetween, composing the potential obstruction
266
267    Returns
268    -------
269    * *bool*
270
271        > True if no obstruction, False else.
272
273    Examples
274    --------
275    >>> import pyvista as pv
276    >>> import pyviewfactor as pvf
277    >>> tri2 = pv.Triangle([[0.0, 1.0, 1.0],[1.0, 1.0, 1.0],[1.0, 0.0, 1.0]])
278    >>> tri1 = pv.Triangle([[1.0, 0.0, 0.0],[1.0, 1.0, 0.0],[0.0, 1.0, 0.0]])
279    >>> obstacle = pv.Circle(radius=3.0)
280    >>> obstacle.translate([0, 0, 0.5], inplace = True)
281    >>> pvf.get_visibility_MT(tri2, tri1, obstacle)
282    False
283
284    """
285
286    # Defining ray origin and vector
287    ray_orig = face1.cell_centers().points[0]
288    ray_end = face2.cell_centers().points[0]
289    ray_dir = ray_end - ray_orig
290    if not obstacle.is_all_triangles:
291        obstacle.triangulate(inplace=True)
292    vis = []
293    for idx in range(obstacle.n_cells):
294        v1, v2, v3 = obstacle.get_cell(idx).points
295        eps = 0.000001
296        edge1 = v2 - v1
297        edge2 = v3 - v1
298        pvec = np.cross(ray_dir, edge2)
299        det = edge1.dot(pvec)
300        if abs(det) < eps:
301            vis.append(True)
302            continue
303        inv_det = 1. / det
304        tvec = ray_orig - v1
305        u = tvec.dot(pvec) * inv_det
306        if u < 0. or u > 1.:
307            vis.append(True)
308            continue
309        qvec = np.cross(tvec, edge1)
310        v = ray_dir.dot(qvec) * inv_det
311        if v < 0. or u + v > 1.:
312            vis.append(True)
313            continue
314        t = edge2.dot(qvec) * inv_det
315        if t < eps:
316            vis.append(True)
317            continue
318        if t > 1.0: # np.dot(ray_dir, ray_dir):
319            vis.append(True)
320            continue
321        vis.append(False)
322        return False
323    return True if all(vis) else False
324
325
326def trunc(values, decs=0):
327    """Return values with *decs* decimals.
328
329    A function to truncate decimals in floats.
330
331    Parameters
332    ----------
333    * **values** : *float*, or *numpy.array* (floats)
334
335        A float value with decimals, or a numpy.array of floats
336
337    * **decs** : *int*, optional
338
339        The number of decimals to keep. The default is 0.
340
341    Returns
342    -------
343    * *float*
344
345        > The same flaot truncated with *decs* decimals, or a the same
346        numpy.array of floats truncated.
347
348    Examples
349    --------
350    >>> import pyvista as pv
351    >>> import pyviewfactor as pvf
352    >>> a = 1.23456789
353    >>> pvf.trunc(a,2)
354    1.23
355    >>> tri1 = pv.Triangle([[0.111111, 1.111111, 1.111111],
356                        [1.222222, 1.222222, 1.222222],
357                        [1.333333, 0.333333, 1.333333]])
358    >>> trunc(tri1.points,2)
359    pyvista_ndarray([[0.11, 1.11, 1.11],
360                     [1.22, 1.22, 1.22],
361                     [1.33, 0.33, 1.33]])
362
363    """
364    return np.trunc(values*10**decs)/(10**decs)
365
366
367@njit  # numba's just in time compilation offers a x2 speedup
368def integrand(x, y, norm_q_carree, norm_p_carree, scal_qpq,
369              scal_qpp, scal_pq, norm_qp_carree):
370    """
371    Return the integrand for a pair of edges of two facets for the view factor
372    computation.
373
374    Used in the *compute_viewfactor* function.
375
376    """
377    integrand_function = np.log(y**2*norm_q_carree
378                                + x**2*norm_p_carree
379                                - 2*y*scal_qpq
380                                + 2*x*scal_qpp
381                                - 2*x*y*scal_pq
382                                + norm_qp_carree
383                                )*scal_pq
384    return integrand_function
385
386
387def compute_viewfactor(cell_1, cell_2, epsilon=0.001, rounding_decimal=6):
388    """
389    View factor computation between cell1 and cell2
390
391    Parameters
392    ----------
393    * **cell_1** : *pyvista.PolyData* facet
394
395        > The first cell.
396
397    * **cell_2** : *pyvista.PolyData* facet
398
399        > The second cell.
400
401    * **epsilon** : *float*
402
403        * Desired precision, default = 0.001
404        * It should *never* be higher than 1e-10, where the precision is close
405          to numpy's intergation error.
406        * A good practice is not to get higher than the number of digits
407          defining your `points`.
408        * A precision  of 1e-8 with points define with 5 digits, for example,
409          will lead to a zero view factor.
410
411    * **rounding_decimal** : *Int*
412
413        > Number of decimals at which rounding shall be done for the points
414          coordinates. This is to avoid numeric issues when the points
415          coordinates are given as for instance 1.234e-13 instead of 0.
416          (default = 6)
417
418    Returns
419    -------
420    * *float*
421
422        > The view factor from **cell_2** to **cell_1**.
423
424    Examples
425    --------
426    >>> import pyvista as pv
427    >>> import pyviewfactor as pvf
428    >>> tri1 = pv.Triangle([[0.0, 1.0, 1.0],[1.0, 1.0, 1.0],[1.0, 0.0, 1.0]])
429    >>> tri2 = pv.Triangle([[1.0, 0.0, 0.0],[1.0, 1.0, 0.0],[0.0, 1.0, 0.0]])
430    >>> pvf.compute_viewfactor(tri1, tri2)
431    0.07665424316999997
432
433    """
434    # View factor initialization
435    FF = 0
436    # Cell 1 preparation
437    cell_1_poly = cell_1
438    cell_1 = cell_1.cast_to_unstructured_grid()
439    # Getting the cell points
440    # cell_1_points = cell_1.get_cell(0).points
441    cell_1_points = np.round(cell_1.get_cell(0).points, rounding_decimal)
442    # Creating vectors describing the oriented edges
443    cell_1_points_roll = np.roll(cell_1_points, -1, axis=0)
444    # [Ai,Bi,Ci,...,Ni] -> [Bi, Ci,..., Ni,Ai]
445    vect_dir_elts_1 = cell_1_points_roll - cell_1_points
446    # [[Bi-Ai], [Ci-Bi],...,[Ni-Ai]]
447    # Ther euse to be a truncature, for test purposes
448    # cell_1_points = trunc(cell_1_points, decs=10)
449
450    # Same for cell 2
451    cell_2_poly = cell_2
452    cell_2 = cell_2.cast_to_unstructured_grid()
453    # cell_2_points = cell_2.get_cell(0).points
454    cell_2_points = np.round(cell_2.get_cell(0).points, rounding_decimal)
455    cell_2_points_roll = np.roll(cell_2_points, -1, axis=0)
456    vect_dir_elts_2 = cell_2_points_roll - cell_2_points
457    #cell_2_points = trunc(cell_2_points, decs=10)
458    # Creation of a matrix N by M, with N the number of points defining cell 1
459    # and M the number of points defining cell 2
460    n_cols = np.shape(cell_2_points)[0]
461    n_rows = np.shape(cell_1_points)[0]
462    # Getting the number of vertexes for eaach cell
463    nb_sommets_1 = n_rows
464    nb_sommets_2 = n_cols
465    # Creation of an empty dataframe of size N by M
466    mat_vectors = np.zeros((n_rows, n_cols))
467    vectors = pd.DataFrame(mat_vectors, dtype=object)
468    # Filling the dataframe
469    for row in range(n_rows):
470        # Getting the coordinates of the vectors starting from Vertex i from
471        # cell 1 to each vertex j from cell 2
472        coord_repeat = np.tile(cell_1_points[row], (nb_sommets_2, 1))
473        vect_sommets_1_to_2 = cell_2_points - coord_repeat
474        # Filling the "vectors" dataframe
475        vectors.iloc[row] = list(vect_sommets_1_to_2)
476    vect_sommets_extra = vectors
477    vect_sommets_intra_1 = vect_dir_elts_1
478    vect_sommets_intra_2 = vect_dir_elts_2
479    # Constants calculations for the contour integral
480    area = cell_2_poly.compute_cell_sizes(area=True)['Area']
481    A_q = area[0]
482    constante = 4*np.pi*A_q
483    # Initialisation of the integration error, and the contribution of each
484    # edge intergal to the overall one
485    err = 0
486    s_i = 0
487    s_j = 0
488    # Getting the number of shared vertexes
489    arr_test = np.argwhere(
490        (cell_2_points[:, None, :] == cell_1_points[:, :]).all(-1)
491    )
492    nbre_sommets_partages = np.shape(arr_test)[0]
493    # If the 2 cells don't share any vertices, "regular" computation of the
494    # integral
495    if nbre_sommets_partages == 0:
496        for n in range(nb_sommets_2):
497            p_n_np1 = tuple(vect_sommets_intra_2[n, :])
498            norm_p_carree = np.dot(p_n_np1, p_n_np1)
499            for m in range(nb_sommets_1):
500                q_m_mp1 = tuple(vect_sommets_intra_1[m, :])
501                norm_q_carree = np.dot(q_m_mp1, q_m_mp1)
502                qm_pn = vect_sommets_extra.loc[m, n]
503                norm_qp_carree = np.dot(qm_pn, qm_pn)
504                scal_qpq = np.dot(qm_pn, q_m_mp1)
505                scal_qpp = np.dot(qm_pn, p_n_np1)
506                scal_pq = np.dot(q_m_mp1, p_n_np1)
507                s_j, err = scipy.integrate.dblquad(
508                    integrand,
509                    0, 1,
510                    lambda x: 0, lambda x: 1,
511                    args=(
512                        norm_q_carree,
513                        norm_p_carree,
514                        scal_qpq,
515                        scal_qpp,
516                        scal_pq,
517                        norm_qp_carree
518                    )
519                )
520                s_i += round(s_j/constante, 11)
521                err += err/(nb_sommets_1 + nb_sommets_2)
522    else:
523        # When the cells share one edge or more:
524        # Apply a virtual displacement by epsilon, according to the
525        # centroid vector
526        centroid_vec = cell_2_poly.cell_centers().points \
527            - cell_1_poly.cell_centers().points
528        # Normalization
529        centroid_vec = centroid_vec / np.sqrt(np.sum(centroid_vec**2))
530        for sommet_j in cell_2_points[:, :]:
531            sommet_j += np.dot(centroid_vec, epsilon)[0]
532        # After the displacement, update the matrix with every vector from
533        # cell 1 vertices to every vertiuces of cell 2
534        for row in range(n_rows):
535            # Getting the coordinates of the vectors starting from Vertex i
536            # from cell 1 to each vertex j from cell 2, as previously
537            coord_repeat = np.tile(cell_1_points[row], (n_cols, 1))
538            vect_sommets_i_to_j = cell_2_points - coord_repeat
539            vectors.iloc[row] = list(vect_sommets_i_to_j)
540        # Then proceed to the regular integral computation
541        for n in range(nb_sommets_2):
542            p_n_np1 = tuple(vect_sommets_intra_2[n, :])
543            norm_p_carree = np.dot(p_n_np1, p_n_np1)
544            for m in range(nb_sommets_1):
545                q_m_mp1 = tuple(vect_sommets_intra_1[m, :])
546                norm_q_carree = np.dot(q_m_mp1, q_m_mp1)
547                qm_pn = vect_sommets_extra.loc[m, n]
548                norm_qp_carree = np.dot(qm_pn, qm_pn)
549                scal_qpq = np.dot(qm_pn, q_m_mp1)
550                scal_qpp = np.dot(qm_pn, p_n_np1)
551                scal_pq = np.dot(q_m_mp1, p_n_np1)
552                s_j, err = scipy.integrate.dblquad(
553                    integrand,
554                    0, 1,
555                    lambda x: 0, lambda x: 1,
556                    args=(
557                        norm_q_carree,
558                        norm_p_carree,
559                        scal_qpq,
560                        scal_qpp,
561                        scal_pq,
562                        norm_qp_carree)
563                )
564                s_i += round(s_j/constante, 11)
565                err += err/(nb_sommets_1 + nb_sommets_2)
566    if s_i > 0:
567        FF = s_i
568    return FF
569
570
571""" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End Of File ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ """