Homec4science
Diffusion CADD_mesher (dev/incorporate-qhull) Browse

rCADDMESH/test/domain_test.cc
58a104d7bd95dev/incorporate-qhull

Last Change
Raw File

domain_test.cc
View Options

#include "cadd_mesher.hh"
#include <iostream>
int main(int argc, char *argv[])
{
const Uint dim = twoD;
// Prepare a lattice
//Real consts[dim] = {1, 1};
Real consts[dim] = {4.04, 4.04};
//Real consts[dim] = {5, 5};
//Real consts[dim] = {2.02, 2.23};
//Real miller[dim*dim] = {1, 1, -1, 1};
Real miller[dim*dim] = {1, 1, -1, 1};
FccLattice<dim> lattice(consts, miller);
std::cout << "Lattice = " << lattice << std::endl;
// Prepare the domain geometries
Real overall_min_max[] = {-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, -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, -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, 19, -29, 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) ; ++i) {
for (int j = 0; j < int(dim) ; ++j) {
refinement_point[0] = (1-2*static_cast<Real>(i))*201;
refinement_point[1] = -static_cast<Real>(j * 301 - (1-j));
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] = -static_cast<Real>(j * 35 - (1-j));
refinement = 0;
domain.setPointRefinement(refinement_point, refinement);
punkte.addPoint(refinement_point);
refinement_point[0] = (1-2*static_cast<Real>(i))*100;
refinement_point[1] = -static_cast<Real>(j * 130 - (1-j));
refinement = 8;
domain.setPointRefinement(refinement_point, refinement);
punkte.addPoint(refinement_point);
aufloesung.push_back(refinement);
}
}
std::cout << punkte << std::endl;
Mesh<dim> testmesh(punkte.getVector(),0);
std::cout << "Text mesh before delaunay = " << std::endl << testmesh;
testmesh.delaunay(10000);
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");
Real point[][2] = {{-30, -5},
{-30, 0},
{-25, 0},
{-25, -5}};
for (Uint i = 0; i < 4 ; ++i) {
Real * pt = point[i];
std::cout <<std::fixed<< "domain.interpolate(" <<pt[0]<< ", "<<pt[1]<< ") = " << domain.interpolate(pt) << std::endl;
}
std::cout << "domain.getPoints() : " << std::endl << domain.getPoints();
try {
Mesh<dim> & main_mesh = domain.getMesh();
main_mesh.dump_paraview("domain_test_output");
main_mesh.dump_msh("domain_test_output_msh");
PointContainer<dim> interface_atoms("interface atoms");
PointContainer<dim> pad_atoms("pad atoms");
Real interface_min_max[] = {-21, 21, -35, 0};
Block<dim> interface(interface_min_max, "refined");
domain.compute_coupling_atoms(interface_atoms, pad_atoms, interface);
interface_atoms.dump_paraview("interface_atoms_2d");
interface_atoms.dump_text("interface_atoms_2d.txt", "contains interface atoms");
pad_atoms.dump_paraview("pad_atoms_2d");
pad_atoms.dump_text("pad_atoms_2d.txt", "contains pad atoms");
std::cout << interface_atoms << std::endl;
domain.dump_atoms_paraview("all_atoms_2d");
main_mesh.computeRadiusRatios();
std::cout << "Mesh quality: q_min = " << main_mesh.getRadiusRatioMin()
<< ", q_max = " << main_mesh.getRadiusRatioMax() << std::endl;
} 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.plot");
return 0;
}
c4science · Help