Skip to content

Commit

Permalink
updates
Browse files Browse the repository at this point in the history
  • Loading branch information
@smartalecH
smartalecH committed Jun 25, 2019
1 parent 6ad454a commit e28c961
Show file tree
Hide file tree
Showing 4 changed files with 37 additions and 51 deletions.
25 changes: 3 additions & 22 deletions python/geom.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -239,28 +239,9 @@ def transform(self, m):
for s in self.H_susceptibilities:
s.transform(m)

def rotate(self, rotations):
for rot in rotations:
if np.count_nonzero(rot) != 1:
raise ValueError("Each rotation vector should only have 1 coordinate.")
if rot.x != 0: # rotate about x axis
self.transform(Matrix(
Vector3(1,0,0),
Vector3(0,np.cos(rot.x),np.sin(rot.x)),
Vector3(0,-np.sin(rot.x),np.cos(rot.x))
))
elif rot.y != 0: # rotate about z axis
self.transform(Matrix(
Vector3(np.cos(rot.y),0,-np.sin(rot.y)),
Vector3(0,1,0),
Vector3(np.sin(rot.y),0,np.cos(rot.y))
))
else:
self.transform(Matrix(
Vector3(np.cos(rot.z),np.sin(rot.z),0),
Vector3(-np.sin(rot.z),np.cos(rot.z),0),
Vector3(0,0,1)
))
def rotate(self, axis, theta):
T = get_rotation_matrix(axis,theta)
self.transform(T)

def epsilon(self,freq):
return self._get_epsmu(self.epsilon_diag, self.epsilon_offdiag, self.E_susceptibilities, self.D_conductivity_diag, self.D_conductivity_offdiag, freq)
Expand Down
10 changes: 6 additions & 4 deletions python/tests/dispersive_eigenmode.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,8 +38,10 @@ def call_chi1(self,material,omega):

def verify_output_and_slice(self,material,omega):
# Since the slice routines average the diagonals, we need to do that too:
chi1 = np.square(np.real(np.sqrt(material.epsilon(omega))))
chi1inv = (np.linalg.inv(chi1))
chi1 = material.epsilon(omega).astype(np.complex128)
if np.any(np.imag(chi1) != 0):
chi1 = np.square(np.real(np.sqrt(chi1)))
chi1inv = np.linalg.inv(chi1)
chi1inv = np.diag(chi1inv)
N = chi1inv.size
n = np.sqrt(N/np.sum(chi1inv))
Expand Down Expand Up @@ -101,7 +103,7 @@ def test_chi1_routine(self):
# Now let's rotate LN
import copy
rotLiNbO3 = copy.deepcopy(LiNbO3)
rotLiNbO3.rotate([mp.Vector3(x=np.radians(28)),mp.Vector3(np.radians(45))])
rotLiNbO3.rotate(mp.Vector3(1,1,1),np.radians(34))
self.call_chi1(rotLiNbO3,w0)
self.call_chi1(rotLiNbO3,w1)

Expand Down Expand Up @@ -138,7 +140,7 @@ def test_get_with_dispersion(self):
# Now let's rotate LN
import copy
rotLiNbO3 = copy.deepcopy(LiNbO3)
rotLiNbO3.rotate([mp.Vector3(x=np.radians(28)),mp.Vector3(np.radians(45))])
rotLiNbO3.rotate(mp.Vector3(1,1,1),np.radians(34))
self.verify_output_and_slice(rotLiNbO3,w0)
self.verify_output_and_slice(rotLiNbO3,w1)

Expand Down
2 changes: 1 addition & 1 deletion src/meep.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ const double nan = -7.0415659787563146e103; // ideally, a value never encountere
class h5file;

// Defined in monitor.cpp
void matrix_invert(std::vector<double> &Vinv, std::vector<double> &V);
void matrix_invert(std::complex<double> (&Vinv)[9], std::complex<double> (&V)[9]);

double pml_quadratic_profile(double, void *);

Expand Down
51 changes: 27 additions & 24 deletions src/monitor.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -225,29 +225,31 @@ double structure::get_chi1inv(component c, direction d, const ivec &origloc, dou
return 0.0;
}

/* Set Vinv = inverse of V, where both V and Vinv are real-symmetric matrices.*/
void matrix_invert(std::vector<double> &Vinv, std::vector<double> &V) {
/* Set Vinv = inverse of V, where both V and Vinv are complex matrices.*/
void matrix_invert(std::complex<double> (&Vinv)[9], std::complex<double> (&V)[9]) {

double det = (V[0 +3*0] * (V[1 + 3*1]*V[2 +3*2] - V[1 + 3*2]*V[2 +3*1]) -
std::complex<double> det = (V[0 +3*0] * (V[1 + 3*1]*V[2 +3*2] - V[1 + 3*2]*V[2 +3*1]) -
V[0 + 3*1] * (V[0 + 3*1]*V[2 + 3*2] - V[1 + 3*2]*V[0 + 3*2]) +
V[0 + 3*2] * (V[0 + 3*1]*V[1 + 3*2] - V[1 + 3*1]*V[0 + 3*2]));

if (det == 0) abort("meep: Matrix is singular, aborting.\n");

Vinv[0 + 3*0] = 1/det * (V[1 + 3*1]*V[2 + 3*2] - V[1 + 3*2]*V[2 + 3*1]);
Vinv[0 + 3*1] = 1/det * (V[0 + 3*2]*V[2 + 3*1] - V[0 + 3*1]*V[2 + 3*2]);
Vinv[0 + 3*2] = 1/det * (V[0 + 3*1]*V[1 + 3*2] - V[0 + 3*2]*V[1 + 3*1]);
Vinv[1 + 3*0] = 1/det * (V[1 + 3*2]*V[2 + 3*0] - V[1 + 3*0]*V[2 + 3*2]);
Vinv[1 + 3*1] = 1/det * (V[0 + 3*0]*V[2 + 3*2] - V[0 + 3*2]*V[2 + 3*0]);
Vinv[1 + 3*2] = 1/det * (V[0 + 3*2]*V[1 + 3*0] - V[0 + 3*0]*V[1 + 3*2]);
Vinv[2 + 3*0] = 1/det * (V[1 + 3*0]*V[2 + 3*1] - V[1 + 3*1]*V[2 + 3*0]);
Vinv[2 + 3*1] = 1/det * (V[0 + 3*1]*V[2 + 3*0] - V[0 + 3*0]*V[2 + 3*1]);
Vinv[2 + 3*2] = 1/det * (V[0 + 3*0]*V[1 + 3*1] - V[0 + 3*1]*V[1 + 3*0]);
if (det.real() == 0 && det.imag() == 0) abort("meep: Matrix is singular, aborting.\n");

Vinv[0 + 3*0] = 1.0/det * (V[1 + 3*1]*V[2 + 3*2] - V[1 + 3*2]*V[2 + 3*1]);
Vinv[0 + 3*1] = 1.0/det * (V[0 + 3*2]*V[2 + 3*1] - V[0 + 3*1]*V[2 + 3*2]);
Vinv[0 + 3*2] = 1.0/det * (V[0 + 3*1]*V[1 + 3*2] - V[0 + 3*2]*V[1 + 3*1]);
Vinv[1 + 3*0] = 1.0/det * (V[1 + 3*2]*V[2 + 3*0] - V[1 + 3*0]*V[2 + 3*2]);
Vinv[1 + 3*1] = 1.0/det * (V[0 + 3*0]*V[2 + 3*2] - V[0 + 3*2]*V[2 + 3*0]);
Vinv[1 + 3*2] = 1.0/det * (V[0 + 3*2]*V[1 + 3*0] - V[0 + 3*0]*V[1 + 3*2]);
Vinv[2 + 3*0] = 1.0/det * (V[1 + 3*0]*V[2 + 3*1] - V[1 + 3*1]*V[2 + 3*0]);
Vinv[2 + 3*1] = 1.0/det * (V[0 + 3*1]*V[2 + 3*0] - V[0 + 3*0]*V[2 + 3*1]);
Vinv[2 + 3*2] = 1.0/det * (V[0 + 3*0]*V[1 + 3*1] - V[0 + 3*1]*V[1 + 3*0]);
}

double structure_chunk::get_chi1inv_at_pt(component c, direction d, int idx, double omega) const {
double res = 0.0;
if (is_mine()){
if (omega == 0)
return chi1inv[c][d] ? chi1inv[c][d][idx] : (d == component_direction(c) ? 1.0 : 0);
// ----------------------------------------------------------------- //
// ---- Step 1: Get instantaneous chi1 tensor ----------------------
// ----------------------------------------------------------------- //
Expand All @@ -268,22 +270,23 @@ double structure_chunk::get_chi1inv_at_pt(component c, direction d, int idx, dou
my_stuff = B_stuff;
}

std::vector<double> chi1_inv_tensor(9,0);
std::vector<double> chi1_tensor(9,0);
std::complex<double> chi1_inv_tensor[9] = {std::complex<double>(1, 0),std::complex<double>(0, 0),std::complex<double>(0, 0),
std::complex<double>(0, 0),std::complex<double>(1, 0),std::complex<double>(0, 0),
std::complex<double>(0, 0),std::complex<double>(0, 0),std::complex<double>(1, 0)
};
std::complex<double> chi1_tensor[9] = {std::complex<double>(1, 0),std::complex<double>(0, 0),std::complex<double>(0, 0),
std::complex<double>(0, 0),std::complex<double>(1, 0),std::complex<double>(0, 0),
std::complex<double>(0, 0),std::complex<double>(0, 0),std::complex<double>(1, 0)
};

// Set up the chi1inv tensor
// Set up the chi1inv tensor with the DC components
for (int com_it=0; com_it<3;com_it++){
for (int dir_int=0;dir_int<3;dir_int++){
if (chi1inv[comp_list[com_it]][dir_int] )
chi1_inv_tensor[com_it + 3*dir_int] = chi1inv[comp_list[com_it]][dir_int][idx];
else if(dir_int == component_direction(comp_list[com_it]))
chi1_inv_tensor[com_it + 3*dir_int] = 1;
else
chi1_inv_tensor[com_it + 3*dir_int] = 0;
}
}


matrix_invert(chi1_tensor, chi1_inv_tensor); // We have the inverse, so let's invert it.

// ----------------------------------------------------------------- //
Expand All @@ -293,7 +296,7 @@ double structure_chunk::get_chi1inv_at_pt(component c, direction d, int idx, dou
// loop over tensor elements
for (int com_it=0; com_it<3;com_it++){
for (int dir_int=0;dir_int<3;dir_int++){
std::complex<double> eps(chi1_tensor[com_it + 3*dir_int],0.0);
std::complex<double> eps = chi1_tensor[com_it + 3*dir_int];
component cc = comp_list[com_it];
direction dd = (direction)dir_int;
// Loop through and add up susceptibility contributions
Expand Down Expand Up @@ -326,7 +329,7 @@ double structure_chunk::get_chi1inv_at_pt(component c, direction d, int idx, dou
// ----------------------------------------------------------------- //

matrix_invert(chi1_inv_tensor, chi1_tensor); // We have the inverse, so let's invert it.
res = chi1_inv_tensor[component_index(c) + 3*d];
res = chi1_inv_tensor[component_index(c) + 3*d].real();
}
return res;
}
Expand Down

0 comments on commit e28c961

Please sign in to comment.