Homec4science
Diffusion CADD_mesher (master) Browse

rCADDMESH/test/domain_test_3d.cc
cfac30bd70a4master

Last Change
Raw File

domain_test_3d.cc
View Options

#include "cadd_mesher.hh"
#include <iostream>
#include <algorithm>
template<typename bla>
void print(bla x) { std::cout << " " << x ;}
int main(int argc, char *argv[])
{
const Uint dim = threeD;
// Prepare a lattice
//Real consts[dim] = {1, 1};
// Real consts[dim] = {4.04, 4.545};
// Real consts[dim] = {5, 5, 5};
Real consts[dim] = {4.04, 4.04, 4.04};
//Real consts[dim] = {2.02, 2.23};
//Real miller[dim*dim] = {1, 1, -1, 1};
Real miller[dim*dim] = {1, 1, 0, -1, 1, 0, 0, 0, 1};
FccLattice<dim> lattice(consts, miller);
// Prepare the domain geometries
Real overall_min_max[] = {-200, 200, -200, 200, -300, 0};
//Real overall_min_max[] = {-4.25, 4.25, -4.25, 4.25};
Block<dim> overall(overall_min_max, "overall");
Real refined_min_max[] = {-25, 25, -25, 25, -35, 0};
// Real refined_min_max[] = {-4.25, 4.25, -4.25, 4.25};
Block<dim> refined(refined_min_max, "refined");
Real atomistic_min_max[] = {-20, 20, -20, 20, -30, 0};
// Real atomistic_min_max[] = {-4.25, 4.25, -4.25, 4.25};
Block<dim> atomistic(atomistic_min_max, "atomistic");
Real atomistic_skinned_min_max[] = {-19.9, 19.9, -19.9, 19.9, -29.9, 0};
// Real atomistic_skinned_min_max[] = {-16, 16, 26, 30};
Block<dim> atomistic_skinned(atomistic_skinned_min_max, "atomistic_skinned");
// Create domain
Domain<dim> domain(overall, refined, atomistic, atomistic_skinned,
lattice, true);
Real refinement_point[dim];
Real refinement;
PointContainer<dim> punkte("FUER DEBUGGING");
std::vector<Real> aufloesung;
// set up auxiliary mesh
for (int i = 0; i < int(dim-1); ++i) {
for (int j = 0; j < int(dim-1); ++j) {
for (int k = 0; k < int(dim-1); ++k) {
refinement_point[0] = (1-2*static_cast<Real>(i))*201;
refinement_point[1] = (1-2*static_cast<Real>(j))*201;
refinement_point[2] = -static_cast<Real>(k * 301 - (1-k));
refinement = 8;
domain.setPointRefinement(refinement_point, refinement);
punkte.addPoint(refinement_point);
aufloesung.push_back(refinement);
refinement_point[0] = (1-2*static_cast<Real>(i))*25;
refinement_point[1] = (1-2*static_cast<Real>(j))*25;
refinement_point[2] = -static_cast<Real>(k * 35 - (1-k));
refinement = 0;
domain.setPointRefinement(refinement_point, refinement);
refinement_point[0] = (1-2*static_cast<Real>(i))*100;
refinement_point[1] = (1-2*static_cast<Real>(j))*100;
refinement_point[2] = -static_cast<Real>(k * 130 - (1-k));
refinement = 4;
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;
domain.getPoints().dump_paraview("points_before_meshing");
domain.getAuxiliaryMesh().dump_paraview("domain_test_auxiliary_mesh_3D");
std::cout << "domain.getPoints() : " << std::endl << domain.getPoints();
try {
Mesh<dim> & main_mesh = domain.getMesh();
main_mesh.dump_paraview("domain_test_output_3D");
main_mesh.dump_msh("domain_test_output_3D.msh");
PointContainer<dim> interface_atoms("interface atoms");
PointContainer<dim> pad_atoms("pad atoms");
Real interface_min_max[] = {-21, 21, -21, 21, -31, 0};
Block<dim> interface(interface_min_max, "refined");
domain.compute_interface_atoms(interface_atoms, pad_atoms, interface);
interface_atoms.dump_paraview("interface_atoms_3d");
interface_atoms.dump_text("interface_atoms_3d.txt", "contains interface atoms");
pad_atoms.dump_paraview("pad_atoms_3d");
pad_atoms.dump_text("pad_atoms_3d.txt", "contains pad atoms");
std::cout << interface_atoms << std::endl;
domain.dump_atoms_paraview("all_atoms_3d");
domain.dump_atoms_lammps("all_atoms_3d.dat");
main_mesh.computeRadiusRatios();
std::cout << "Mesh quality: q_min = " << main_mesh.getRadiusRatioMin()
<< ", q_max = " << main_mesh.getRadiusRatioMax() << std::endl;
std::vector<int> degenerates;
main_mesh.getDegenerateElements(10, degenerates);
std::cout << std::endl;
for_each(degenerates.begin(), degenerates.end(), print<int>);
std::cout << std::endl;
std::cout << "quality of element 1325: " << main_mesh.computeRadiusRatios(degenerates[0])[0] << std::endl;
domain.splitDegenerates(10);
Mesh<dim> & new_main_mesh = domain.getMesh();
new_main_mesh.computeRadiusRatios();
std::cout << "Mesh quality after split: q_min = " << new_main_mesh.getRadiusRatioMin()
<< ", q_max = " << new_main_mesh.getRadiusRatioMax() << std::endl;
new_main_mesh.dump_paraview("domain_test_output_3D_after_split");
//new_main_mesh.dump_msh("domain_test_output_msh_3D_after_split");
} catch (std::string & error) {
std::cout << "There's been an error during meshing. You can consult "
<< "'tree_diagnostic' for debugging" << std::endl;
std::cout << error << std::endl;
}
domain.getTree().dump("tree_3D.plot");
return 0;
}
c4science · Help