# Path: openmodes/operator/singularities.py
# def singular_impedance_rwg(basis, num_terms, rel_tol, normals):
#     """Precalculate the singular impedance terms for an object
# 
#     Parameters
#     ----------
#     basis: LinearTriangleBasis object
#         The basis functions representing the object for which to calculate the
#         singularities
#     num_terms: integer
#         The number of singular terms to extract
#     rel_tol: float
#         The desired relative tolerance of the singular integrals
#     normals: ndarray(num_triangles, 3) of float, optional
#         The surface normals, required for n x operator forms
# 
#     Returns
#     -------
#     singular_terms : dictionary of SingularSparse object
#         Each key "T_EFIE", "T_MFIE" and "N_MFIE" has a value withthe sparse
#         array of singular impedance terms for that operator
#     """
# 
#     if not isinstance(basis, LinearTriangleBasis):
#         raise ValueError("Basis functions are not RWG based")
# 
#     # Check if this part's singularities have previously been calculated
#     # Note that higher accuracy calculations will not be used if a lower
#     # accuracy is requested. This avoids non-deterministic behaviour.
#     normals_hash = hashlib.sha1(normals).hexdigest()
#     unique_id = ("RWG", basis.id, rel_tol, normals_hash, num_terms)
#     if unique_id in cached_singular_terms:
#         return cached_singular_terms[unique_id]
# 
#     sharing_nodes = basis.mesh.triangles_sharing_nodes()
# 
#     # slightly inefficient reordering and resizing of nodes array
#     polygons = np.ascontiguousarray(basis.mesh.polygons)
#     nodes = np.ascontiguousarray(basis.mesh.nodes, dtype=np.float64)
#     num_faces = len(polygons)
# 
#     singular_terms = {"T_EFIE": MultiSparse([(np.float64, (num_terms,)),  # phi
#                                              (np.float64, (num_terms, 3, 3))]),  # A
#                       "T_MFIE": MultiSparse([(np.float64, (num_terms, 3, 3))]),  # A
#                       "N_MFIE": MultiSparse([(np.float64, (num_terms, 3, 3))]),  # A
#                       }
# 
#     logging.info("Integrating singular terms for basis function %s, with %d "
#                  "terms, relative tolerance %e" % (basis, num_terms, rel_tol))
# 
#     # find the neighbouring triangles (including self terms) to integrate
#     # singular part
#     for p in range(0, num_faces):  # observer
#         sharing_triangles = set()
#         for node in polygons[p]:
#             sharing_triangles = sharing_triangles.union(sharing_nodes[node])
#         # find any neighbouring elements which are touching
#         for q in sharing_triangles:  # source
#             # at least one node is shared between p and q
#             normal = normals[p]
#             res = face_integrals_yla_oijala(nodes[polygons[q]], rule.points, rule.weights,
#                                             nodes[polygons[p]], normal_o=normal)
#             if q != p:
#                 # The self triangle terms are not evaluated for MFIE
#                 singular_terms["T_MFIE"][p, q] = (res[3],)
#                 singular_terms["N_MFIE"][p, q] = (res[2],)
#             singular_terms["T_EFIE"][p, q] = (res[1], res[0])
# 
#     # Arrays are currently put into fortran order, under the assumption
#     # that they will mostly be used by fortran routines.
#     cached_singular_terms[unique_id] = {k: v.to_csr(order='F')
#                                         for k, v in singular_terms.items()}
#     return cached_singular_terms[unique_id]
# 
# Path: openmodes/constants.py



import numpy as np
from openmodes.operator.singularities import singular_impedance_rwg
from openmodes.core import z_mfie_faces_self, z_mfie_faces_mutual
from openmodes.core import z_efie_faces_self, z_efie_faces_mutual
from openmodes.constants import pi, c
                                integration_rule.weights, *singular_terms)

        transform_L_s = transform_L_o
        transform_S_s = transform_S_o

    else:
        # calculate mutual impedance

        num_faces_s = len(basis_s.mesh.polygons)

        res = z_efie_faces_mutual(nodes_o, basis_o.mesh.polygons, nodes_s,
                                  basis_s.mesh.polygons, gamma_0,
                                  integration_rule.points,
                                  integration_rule.weights)

        transform_L_s, transform_S_s = basis_s.transformation_matrices

    A_faces, phi_faces, A_dgamma_faces, phi_dgamma_faces = res

    if np.any(np.isnan(A_faces)) or np.any(np.isnan(phi_faces)):
        raise ValueError("NaN returned in impedance matrix")

    if np.any(np.isnan(A_dgamma_faces)) or np.any(np.isnan(phi_dgamma_faces)):
        raise ValueError("NaN returned in impedance matrix derivative")

    L = transform_L_o.dot(transform_L_s.dot(A_faces.reshape(num_faces_o*3,
                                                            num_faces_s*3,
                                                            order='C').T).T)
    S = transform_S_o.dot(transform_S_s.dot(phi_faces.T).T)


# Path: openmodes/operator/singularities.py
# def singular_impedance_rwg(basis, num_terms, rel_tol, normals):
#     """Precalculate the singular impedance terms for an object
# 
#     Parameters
#     ----------
#     basis: LinearTriangleBasis object
#         The basis functions representing the object for which to calculate the
#         singularities
#     num_terms: integer
#         The number of singular terms to extract
#     rel_tol: float
#         The desired relative tolerance of the singular integrals
#     normals: ndarray(num_triangles, 3) of float, optional
#         The surface normals, required for n x operator forms
# 
#     Returns
#     -------
#     singular_terms : dictionary of SingularSparse object
#         Each key "T_EFIE", "T_MFIE" and "N_MFIE" has a value withthe sparse
#         array of singular impedance terms for that operator
#     """
# 
#     if not isinstance(basis, LinearTriangleBasis):
#         raise ValueError("Basis functions are not RWG based")
# 
#     # Check if this part's singularities have previously been calculated
#     # Note that higher accuracy calculations will not be used if a lower
#     # accuracy is requested. This avoids non-deterministic behaviour.
#     normals_hash = hashlib.sha1(normals).hexdigest()
#     unique_id = ("RWG", basis.id, rel_tol, normals_hash, num_terms)
#     if unique_id in cached_singular_terms:
#         return cached_singular_terms[unique_id]
# 
#     sharing_nodes = basis.mesh.triangles_sharing_nodes()
# 
#     # slightly inefficient reordering and resizing of nodes array
#     polygons = np.ascontiguousarray(basis.mesh.polygons)
#     nodes = np.ascontiguousarray(basis.mesh.nodes, dtype=np.float64)
#     num_faces = len(polygons)
# 
#     singular_terms = {"T_EFIE": MultiSparse([(np.float64, (num_terms,)),  # phi
#                                              (np.float64, (num_terms, 3, 3))]),  # A
#                       "T_MFIE": MultiSparse([(np.float64, (num_terms, 3, 3))]),  # A
#                       "N_MFIE": MultiSparse([(np.float64, (num_terms, 3, 3))]),  # A
#                       }
# 
#     logging.info("Integrating singular terms for basis function %s, with %d "
#                  "terms, relative tolerance %e" % (basis, num_terms, rel_tol))
# 
#     # find the neighbouring triangles (including self terms) to integrate
#     # singular part
#     for p in range(0, num_faces):  # observer
#         sharing_triangles = set()
#         for node in polygons[p]:
#             sharing_triangles = sharing_triangles.union(sharing_nodes[node])
#         # find any neighbouring elements which are touching
#         for q in sharing_triangles:  # source
#             # at least one node is shared between p and q
#             normal = normals[p]
#             res = face_integrals_yla_oijala(nodes[polygons[q]], rule.points, rule.weights,
#                                             nodes[polygons[p]], normal_o=normal)
#             if q != p:
#                 # The self triangle terms are not evaluated for MFIE
#                 singular_terms["T_MFIE"][p, q] = (res[3],)
#                 singular_terms["N_MFIE"][p, q] = (res[2],)
#             singular_terms["T_EFIE"][p, q] = (res[1], res[0])
# 
#     # Arrays are currently put into fortran order, under the assumption
#     # that they will mostly be used by fortran routines.
#     cached_singular_terms[unique_id] = {k: v.to_csr(order='F')
#                                         for k, v in singular_terms.items()}
#     return cached_singular_terms[unique_id]
# 
# Path: openmodes/constants.py



# Path: openmodes/operator/rwg.py
import numpy as np
from openmodes.operator.singularities import singular_impedance_rwg
from openmodes.core import z_mfie_faces_self, z_mfie_faces_mutual
from openmodes.core import z_efie_faces_self, z_efie_faces_mutual
from openmodes.constants import pi, c

                                integration_rule.weights, *singular_terms)

        transform_L_s = transform_L_o
        transform_S_s = transform_S_o

    else:
        # calculate mutual impedance

        num_faces_s = len(basis_s.mesh.polygons)

        res = z_efie_faces_mutual(nodes_o, basis_o.mesh.polygons, nodes_s,
                                  basis_s.mesh.polygons, gamma_0,
                                  integration_rule.points,
                                  integration_rule.weights)

        transform_L_s, transform_S_s = basis_s.transformation_matrices

    A_faces, phi_faces, A_dgamma_faces, phi_dgamma_faces = res

    if np.any(np.isnan(A_faces)) or np.any(np.isnan(phi_faces)):
        raise ValueError("NaN returned in impedance matrix")

    if np.any(np.isnan(A_dgamma_faces)) or np.any(np.isnan(phi_dgamma_faces)):
        raise ValueError("NaN returned in impedance matrix derivative")

    L = transform_L_o.dot(transform_L_s.dot(A_faces.reshape(num_faces_o*3,
                                                            num_faces_s*3,
                                                            order='C').T).T)
    S = transform_S_o.dot(transform_S_s.dot(phi_faces.T).T)


    L /= 4*pi