Homec4science
Diffusion CADD_mesher (master) Browse

rCADDMESH/test/rotation_test_3d.cc
58a104d7bd95master

Last Change
Raw File

rotation_test_3d.cc
View Options

#include "cadd_mesher.hh"
#include <iostream>
#include <algorithm>
#include <cmath>
template<typename bla>
void print(bla x) { std::cout << " " << x ;}
int main(int argc, char *argv[])
{
const Uint dim = threeD;
// Prepare a lattice
Real a = 4.04;
Real consts[dim] = {a, a, a};
Real miller[dim*dim] = {-1, 1, 0, 1, 1, 1, 1, 1, -2};
FccLattice<dim> lattice(consts, miller);
a *= 2;
Real x_min = 0, y_min = 0, z_min = 0;
Real x_max = std::hypot(miller[0], std::hypot(miller[1], miller[2])) * a;
Real y_max = std::hypot(miller[3], std::hypot(miller[4], miller[5])) * a;
Real z_max = std::hypot(miller[6], std::hypot(miller[7], miller[8])) * a;
// max coord
Real coord_max = 0;
coord_max = std::fmax(coord_max, std::fabs(x_min));
coord_max = std::fmax(coord_max, std::fabs(y_min));
coord_max = std::fmax(coord_max, std::fabs(z_min));
coord_max = std::fmax(coord_max, std::fabs(x_max));
coord_max = std::fmax(coord_max, std::fabs(y_max));
coord_max = std::fmax(coord_max, std::fabs(z_max));
// Prepare the domain geometries
Real overall_min_max[] = {x_min, x_max, y_min, y_max, z_min, z_max};
Block<dim> overall(overall_min_max, "overall");
// Create domain
bool boundary_special_treatment = false;
Domain<dim> domain(overall, overall, overall, overall,
lattice, boundary_special_treatment);
Real refinement_point[dim];
Real refinement;
PointContainer<dim> punkte("FUER DEBUGGING");
std::vector<Real> aufloesung;
// set up auxiliary mesh
for (int i = -1; i < 2; i+=2) {
refinement_point[0] = i * coord_max;
for (int j = -1; j < 2; j+=2) {
refinement_point[1] = j * coord_max;
for (int k = -1; k < 2; k+=2) {
refinement_point[2] = k * coord_max;
refinement = 0;
domain.setPointRefinement(refinement_point, refinement);
punkte.addPoint(refinement_point);
aufloesung.push_back(refinement);
std::cout << "no of points = " << punkte.getSize() << std::endl;
}
}
}
std::cout << punkte << std::endl;
// Mesh<dim> testmesh(punkte.getVector());
// testmesh.delaunay(1000);
// testmesh.initInterpolation();
// std::cout << "Test mesh = " << std::endl << testmesh;
// do all the domain test stuff here
try {
domain.fill_points();
} catch (std::string & error) {
std::cout << error << std::endl;
}
std::cout << "number of tree boxes = " << domain.getTree().getNbBoxes() << std:: endl;
PointContainer<dim> atoms = domain.getAtoms();
std::vector<Real> bounds = atoms.computeBounds();
Block<dim> atom_geom(bounds, "atomic geom");
Real offset[dim] = {a/10, a/10, a/10};
atom_geom.shift(offset);
Real tol = 1e-4;
atoms.filter(atom_geom, tol);
atoms.dump_paraview("rot_atoms");
atoms.dump_lammps("rot_atoms.dat");
return 0;
}
c4science · Help