Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88765466
lattice.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Oct 20, 14:07
Size
6 KB
Mime Type
text/x-c++
Expires
Tue, Oct 22, 14:07 (1 d, 16 h)
Engine
blob
Format
Raw Data
Handle
21778173
Attached To
rCADDMESH CADD_mesher
lattice.cc
View Options
/**
* @file lattice.cc
* @author Till Junge <junge@lsmspc42.epfl.ch>
* @date Fri Apr 13 15:55:56 2012
*
* @brief
*
* @section LICENSE
*
* <insert lisence here>
*
*/
/* -------------------------------------------------------------------------- */
#include "lattice.hh"
#include <gsl/gsl_blas.h>
#include "gsl/gsl_vector.h"
#include "gsl/gsl_matrix.h"
#include <string>
#include <limits>
#include <sstream>
#include <cmath>
template<Uint DIM>
Real determinant(gsl_matrix * rotation) throw(std::string){
Real determinant = 0;
Uint num_terms;
if (DIM == 2) {
num_terms = 1;
} else if (DIM == 3) {
num_terms = 3;
} else {
std::stringstream error_stream ;
error_stream << "determinant for matrices other than 2x2 or 3x3 are not "
<< "implemented" << std::endl;
std::cout << error_stream.str();
throw(error_stream.str());
}
for (Uint i = 0; i < num_terms ; ++i) {
Real increment = 1;
Real decrement = 1;
for (Uint j = 0; j < DIM ; ++j) {
increment *= gsl_matrix_get(rotation, j, (j+i)%DIM);
decrement *= gsl_matrix_get(rotation, DIM-1-j, (j+i)%DIM);
}
determinant += increment;
determinant -= decrement;
}
return determinant;
}
template<Uint DIM>
Lattice<DIM>::Lattice(Real * _constants, Real * miller) throw(std::string)
:atoms("atoms")
{
this->constants[DIM] = 1;
for (Uint i = 0; i < DIM; ++i) {
this->constants[i] = _constants[i];
this->constants[DIM] *= _constants[i];
}
this->constants[DIM] = pow(this->constants[DIM], 1/static_cast<Real>(DIM));
this->millers = gsl_matrix_calloc(DIM, DIM);
this->rotation_matrix = gsl_matrix_calloc(DIM, DIM);
if (miller == NULL) {
gsl_vector_view diag = gsl_matrix_diagonal(this->millers);
gsl_vector_set_all(&diag.vector, 1);
diag = gsl_matrix_diagonal(this->rotation_matrix);
gsl_vector_set_all(&diag.vector, 1);
} else {
for (Uint i = 0; i < DIM; ++i) {
for (Uint j = 0; j < DIM; ++j) {
gsl_matrix_set(this->millers, i, j, miller[i*DIM + j]);
}
gsl_vector_const_view miller_row =
gsl_matrix_const_row(this->millers, i);
Real norm = gsl_blas_dnrm2(&miller_row.vector);
gsl_vector_view rotation_column =
gsl_matrix_column(this->rotation_matrix, i);
gsl_vector_memcpy(&rotation_column.vector, &miller_row.vector);
gsl_vector_scale(&rotation_column.vector, 1./norm);
}
}
Real det = determinant<DIM>(this->rotation_matrix);
if (fabs(det - 1) > std::numeric_limits<Real>::epsilon()){
std::stringstream error_stream;
error_stream << "Error: the miller indices you specified yield an invalid "
<< "rotation matrix (The determinant is " << det << " instead "
<< "of 1)." << std::endl;
for (Uint i = 0; i < DIM ; ++i) {
error_stream << "[";
for (Uint j = 0; j < DIM-1 ; ++j) {
error_stream << gsl_matrix_get(this->rotation_matrix, i, j) << ", ";
}
error_stream << gsl_matrix_get(this->rotation_matrix, i, DIM-1) << "]"
<< std::endl;
}
error_stream <<" Make sure they for a normal basis and comply with "
<< "the right hand rule. Lattice details: "
<< std::endl;
this->printself(error_stream, indent_incr);
std::cout << error_stream.str();
throw error_stream.str();
}
}
template<Uint DIM>
Lattice<DIM>::Lattice(const Lattice<DIM> & other) {
for (Uint i = 0 ; i < DIM+1 ; ++i) {
this->constants[i] = other.constants[i];
}
this->millers = gsl_matrix_calloc(DIM, DIM);
gsl_matrix_memcpy (this->millers, other.millers);
this->rotation_matrix = gsl_matrix_calloc(DIM, DIM);
gsl_matrix_memcpy(this->rotation_matrix, other.rotation_matrix);
this->atoms = other.atoms;
}
template<Uint DIM>
Lattice<DIM>::~Lattice<DIM>()
{
gsl_matrix_free(this->millers);
gsl_matrix_free(this->rotation_matrix);
}
template<Uint DIM>
void Lattice<DIM>::fillAtoms(int * lattice_coords, PointContainer<DIM> & container)
{
gsl_vector * current_point = gsl_vector_alloc(DIM);
gsl_vector * rotated_point = gsl_vector_alloc(DIM);
for (Uint i = 0; i < this->atoms.getSize(); ++i) {
PointRef<DIM> current_atom = this->atoms.getPoint(i);
for (Uint j = 0; j < DIM; ++j) {
gsl_vector_set(current_point, j, current_atom.getComponent(j) +
lattice_coords[j] * this->constants[j]);
}
gsl_blas_dgemv (CblasNoTrans,
1.0, this->rotation_matrix, current_point,
0.0, rotated_point);
container.addPoint(rotated_point->data);
}
gsl_vector_free(current_point);
gsl_vector_free(rotated_point);
}
template<Uint DIM>
void Lattice<DIM>::fillAtoms(Real * sw_coords, PointContainer<DIM> & container)
{
gsl_vector * current_point = gsl_vector_alloc(DIM);
gsl_vector * rotated_point = gsl_vector_alloc(DIM);
for (Uint i = 0; i < this->atoms.getSize(); ++i) {
PointRef<DIM> current_atom = this->atoms.getPoint(i);
for (Uint j = 0; j < DIM; ++j) {
gsl_vector_set(current_point, j, current_atom.getComponent(j));
}
gsl_blas_dgemv (CblasNoTrans,
1.0, this->rotation_matrix, current_point,
0.0, rotated_point);
for (Uint i = 0; i < DIM ; ++i) {
rotated_point->data[i] += sw_coords[i];
}
container.addPoint(rotated_point->data);
}
gsl_vector_free(current_point);
gsl_vector_free(rotated_point);
}
template <Uint DIM>
void Lattice<DIM>::printself(std::ostream & stream, int indent) const
{
indent_space(stream, indent);
stream << this->lattice_type << DIM << "D Lattice, a = (";
for (Uint i = 0; i < DIM-1; ++i) {
stream << this->constants[i] << ", ";
}
stream << this->constants[DIM-1] << "):" << std::endl;
this->atoms.printself(stream, indent+indent_incr);
indent_space(stream, indent+indent_incr);
stream << "Miller orientations :" << std::endl;
for (Uint i = 0; i < DIM; ++i) {
indent_space(stream, indent + 2*indent_incr);
stream << "[";
for (Uint j = 0; j < DIM-1; ++j) {
stream << gsl_matrix_get(this->millers, i, j) << ", ";
}
stream << gsl_matrix_get(this->millers, i, DIM-1) << "]" << std::endl;
}
}
template <Uint DIM>
std::string Lattice<DIM>::lattice_type = "";
template class Lattice<2>;
template class Lattice<3>;
Event Timeline
Log In to Comment