# Path: sucemfem/Utilities/MeshGenerators.py
# def get_centred_cube(domain_size, max_edge_len, centred_element_coordinate=None):
#     """ Generate a cube mesh centred arount [0,0,0]
# 
#     optionally translates mesh slightly such that
#     centred_element_coordinate is at the centre of an element.
#     """
#     domain_subdivisions = N.array(N.ceil(N.sqrt(2)*domain_size/max_edge_len), N.uint)
#     mesh = dolfin.UnitCube(*domain_subdivisions)
#     # Transform mesh to correct dimensions
#     mesh.coordinates()[:] *= domain_size
#     mesh.coordinates()[:] -= domain_size/2
#     if centred_element_coordinate is not None:
#         ## Translate mesh slightly so that source coordinate lies at
#         ## centroid of an element
#         centred_element_point = dolfin.Point(*centred_element_coordinate)
# #        source_elnos = mesh.all_intersected_entities(centred_element_point)
#         source_elnos = mesh.intersected_cells(centred_element_point)
#         closest_elno = source_elnos[(N.argmin(
#             [centred_element_point.distance(dolfin.Cell(mesh, i).midpoint())
#              for i in source_elnos]))]
#         centre_pt = dolfin.Cell(mesh, closest_elno).midpoint()
#         centre_coord = N.array([centre_pt.x(), centre_pt.y(), centre_pt.z()])
#         # The mesh intersect operator caches certain data
#         # structures. We need to clear them if the mesh coordinates
#         # are changed after calling any mesh intersection methods.
#         mesh.intersection_operator().clear()
#         mesh.coordinates()[:] -= centre_coord
#     return mesh
# 
# Path: sucemfem/Consts.py
# Z0 = mu0*c0                              
# 
# Path: sucemfem/Sources/fillament_current_source.py
# class FillamentCurrentSource(CurrentSource):
#     no_integration_points = 100
# 
#     def __init__(self, *names, **kwargs):
#         self._dirty = True
#         
#     def set_no_integration_points(self, no_integration_points):
#         """Set number of integration points to use along the length of
#         the fillament"""
#         self.no_integration_points = no_integration_points
#     
#     def set_source_endpoints(self, source_endpoints):
#         """Set the current filament endpoints.
# 
#         @param source_endpoints: 2x3 array with the coordinates of the start
#             and end point of the current fillament. The conventional
#             current flows from the start point towards the endpoint.
#         """
#         self.source_endpoints = source_endpoints
#         self._dirty = True
# 
#     def set_value(self, value):
#         """Set line current value in Amperes"""
#         self.value = value
#         self._dirty = True
#         
#     def _update(self):
#         if self._dirty:
#             self.source_start, self.source_end = self.source_endpoints
#             self.source_delta = self.source_end - self.source_start
#             self.vector_value = unit_vector(self.source_delta)*self.value
#         self._dirty == False
#             
#     def get_contribution(self):
#         self._update()
#         source_len = vector_length(self.source_delta)
#         no_pts = self.no_integration_points
#         if no_pts == 1:
#             intg_pts = [self.source_start + self.source_delta*0.5]
#         else:
#             intg_pts = self.source_start + self.source_delta*np.linspace(
#                 0,1,no_pts)[:, np.newaxis]
# 
#         contribs = collections.defaultdict(lambda : 0.)
#         point_magnitude = self.vector_value*source_len/no_pts
#         for pt in intg_pts:
#             dnos, vals = calc_pointsource_contrib(
#                 self.function_space, pt, point_magnitude)
#             for dn, v in zip(dnos, vals):
#                 contribs[dn] = contribs[dn] + v
# 
#         dofnos = np.array(contribs.keys(), dtype=np.uint)
#         rhs_contribs = np.array(contribs.values())
#         return dofnos, rhs_contribs


import dolfin
import numpy as np
import pickle
import sys
import sucemfem.Utilities.Optimization
from sucemfem.Utilities.MeshGenerators import get_centred_cube
from sucemfem.Consts import c0
from sucemfem.Sources.fillament_current_source import FillamentCurrentSource



sucemfem.Utilities.Optimization.set_dolfin_optimisation(True)


## Postprocessing requests
theta_deg = np.linspace(0, 180, 181)
no_ff_pts = len(theta_deg)
phi_deg = np.zeros(no_ff_pts)


## Problem parameters
freq =  1.0e+9                          # Frequency
lam = c0/freq
l = lam/4                               # Dipole length
# l = 2.3*lam                               # Dipole length
I = 1.0                                 # Dipole current
source_direction = np.array([0,0,1.])    # Source orientation
source_centre = np.array([0,0,0.])        # Position of the source
source_endpoints =  np.array(
    [-source_direction*l/2, source_direction*l/2]) + source_centre

## Discretisation settings
order = 2
domain_size = np.array([lam]*3)/2
max_edge_len = lam/3
mesh = get_centred_cube(domain_size, max_edge_len)
## Implementation
V = dolfin.FunctionSpace(mesh, "Nedelec 1st kind H(curl)", order)

# Path: sucemfem/Utilities/MeshGenerators.py
# def get_centred_cube(domain_size, max_edge_len, centred_element_coordinate=None):
#     """ Generate a cube mesh centred arount [0,0,0]
# 
#     optionally translates mesh slightly such that
#     centred_element_coordinate is at the centre of an element.
#     """
#     domain_subdivisions = N.array(N.ceil(N.sqrt(2)*domain_size/max_edge_len), N.uint)
#     mesh = dolfin.UnitCube(*domain_subdivisions)
#     # Transform mesh to correct dimensions
#     mesh.coordinates()[:] *= domain_size
#     mesh.coordinates()[:] -= domain_size/2
#     if centred_element_coordinate is not None:
#         ## Translate mesh slightly so that source coordinate lies at
#         ## centroid of an element
#         centred_element_point = dolfin.Point(*centred_element_coordinate)
# #        source_elnos = mesh.all_intersected_entities(centred_element_point)
#         source_elnos = mesh.intersected_cells(centred_element_point)
#         closest_elno = source_elnos[(N.argmin(
#             [centred_element_point.distance(dolfin.Cell(mesh, i).midpoint())
#              for i in source_elnos]))]
#         centre_pt = dolfin.Cell(mesh, closest_elno).midpoint()
#         centre_coord = N.array([centre_pt.x(), centre_pt.y(), centre_pt.z()])
#         # The mesh intersect operator caches certain data
#         # structures. We need to clear them if the mesh coordinates
#         # are changed after calling any mesh intersection methods.
#         mesh.intersection_operator().clear()
#         mesh.coordinates()[:] -= centre_coord
#     return mesh
# 
# Path: sucemfem/Consts.py
# Z0 = mu0*c0                              
# 
# Path: sucemfem/Sources/fillament_current_source.py
# class FillamentCurrentSource(CurrentSource):
#     no_integration_points = 100
# 
#     def __init__(self, *names, **kwargs):
#         self._dirty = True
#         
#     def set_no_integration_points(self, no_integration_points):
#         """Set number of integration points to use along the length of
#         the fillament"""
#         self.no_integration_points = no_integration_points
#     
#     def set_source_endpoints(self, source_endpoints):
#         """Set the current filament endpoints.
# 
#         @param source_endpoints: 2x3 array with the coordinates of the start
#             and end point of the current fillament. The conventional
#             current flows from the start point towards the endpoint.
#         """
#         self.source_endpoints = source_endpoints
#         self._dirty = True
# 
#     def set_value(self, value):
#         """Set line current value in Amperes"""
#         self.value = value
#         self._dirty = True
#         
#     def _update(self):
#         if self._dirty:
#             self.source_start, self.source_end = self.source_endpoints
#             self.source_delta = self.source_end - self.source_start
#             self.vector_value = unit_vector(self.source_delta)*self.value
#         self._dirty == False
#             
#     def get_contribution(self):
#         self._update()
#         source_len = vector_length(self.source_delta)
#         no_pts = self.no_integration_points
#         if no_pts == 1:
#             intg_pts = [self.source_start + self.source_delta*0.5]
#         else:
#             intg_pts = self.source_start + self.source_delta*np.linspace(
#                 0,1,no_pts)[:, np.newaxis]
# 
#         contribs = collections.defaultdict(lambda : 0.)
#         point_magnitude = self.vector_value*source_len/no_pts
#         for pt in intg_pts:
#             dnos, vals = calc_pointsource_contrib(
#                 self.function_space, pt, point_magnitude)
#             for dn, v in zip(dnos, vals):
#                 contribs[dn] = contribs[dn] + v
# 
#         dofnos = np.array(contribs.keys(), dtype=np.uint)
#         rhs_contribs = np.array(contribs.values())
#         return dofnos, rhs_contribs


# Path: sucemfem/Sources/tests/generate_fillament_source_data.py
import dolfin
import numpy as np
import pickle
import sys
import sucemfem.Utilities.Optimization
from sucemfem.Utilities.MeshGenerators import get_centred_cube
from sucemfem.Consts import c0
from sucemfem.Sources.fillament_current_source import FillamentCurrentSource




sucemfem.Utilities.Optimization.set_dolfin_optimisation(True)


## Postprocessing requests
theta_deg = np.linspace(0, 180, 181)
no_ff_pts = len(theta_deg)
phi_deg = np.zeros(no_ff_pts)


## Problem parameters
freq =  1.0e+9                          # Frequency
lam = c0/freq
l = lam/4                               # Dipole length
# l = 2.3*lam                               # Dipole length
I = 1.0                                 # Dipole current
source_direction = np.array([0,0,1.])    # Source orientation
source_centre = np.array([0,0,0.])        # Position of the source
source_endpoints =  np.array(
    [-source_direction*l/2, source_direction*l/2]) + source_centre

## Discretisation settings
order = 2
domain_size = np.array([lam]*3)/2
max_edge_len = lam/3
mesh = get_centred_cube(domain_size, max_edge_len)
## Implementation
V = dolfin.FunctionSpace(mesh, "Nedelec 1st kind H(curl)", order)

dipole_source = FillamentCurrentSource()