Homec4science
Diffusion CADD_mesher (master) Browse

rCADDMESH/src/domain/mini_mesh.cc
9102204a367dmaster

Last Change
Raw File

mini_mesh.cc
View Options

#include "mini_mesh.hh"
#include "mesh_helper.hh"
#include "common.hh"
#include <cmath>
/* -------------------------------------------------------------------------- */
template <Uint DIM>
MiniMesh<DIM>::MiniMesh(const std::vector<const int*> & connectivity,
const std::vector<int> & free_nodes,
PointContainer<DIM> & nodes)
:connectivity(connectivity), free_nodes(free_nodes), nodes(nodes),
grad(DIM*free_nodes.size()) {}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
Real MiniMesh<DIM>::computeQuality() {
Real quality = 1;
for (auto elem: this->connectivity) {
keep_max(quality,
computeRadiusRatio<DIM>(this->nodes.getArray(), elem,
this->nodes.getSize()));
}
return quality;
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void MiniMesh<DIM>::updateGradient(const Real & step_size) {
const Real ref_quality = this->computeQuality();
Uint nb_free_nodes = this->free_nodes.size();
for (Uint local_nd_id = 0 ; local_nd_id < nb_free_nodes ; ++local_nd_id) {
Uint node_id = this->free_nodes[local_nd_id];
for (Uint i = 0 ; i < DIM ; ++i) {
this->nodes[node_id][i] += step_size;
const Real new_quality = this->computeQuality();
this->nodes[node_id][i] -= step_size;
this->grad.at(DIM*local_nd_id+i) = (new_quality-ref_quality)/step_size;
}
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void MiniMesh<DIM>::normGradient(Real step_size) {
Real max_step = 0;
Uint nb_free_nodes = this->free_nodes.size();
for (Uint step_id = 0 ; step_id < nb_free_nodes ; ++step_id) {
keep_max(max_step, vector_norm<DIM>(&this->grad[DIM*step_id]));
}
Real factor = step_size/max_step;
for (Real & val: this->grad) {
val *= factor;
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
Uint MiniMesh<DIM>::relax(Real initial_step_size,
Uint max_iter, Real threshold,
Real min_step_size,
bool verbose) {
Uint iter = 0;
Real quality = this->computeQuality();
Real old_quality = quality + 1;
Uint nb_free_nodes = this->free_nodes.size();
Uint active_max_iter = max_iter;
while (quality > threshold) {
if (verbose) {
std::cout << " At beginning of iteration " << iter+1
<< ", local quality is " << quality << std::endl;
}
if(++iter >= active_max_iter) {
throw StopIterError("reached max_iter without reaching threshold");
}
this->updateGradient(initial_step_size);
this->normGradient(initial_step_size);
for (Uint local_nd_id = 0 ; local_nd_id < nb_free_nodes ; ++local_nd_id) {
Uint global_nd_id = this->free_nodes[local_nd_id];
for (Uint i = 0 ; i < DIM ; ++i) {
this->nodes[global_nd_id][i] -= this->grad[DIM*local_nd_id+i];
}
}
old_quality = quality;
quality = this->computeQuality();
// make sure it was an improvement
if ((old_quality < quality) && (initial_step_size > min_step_size)) {
quality = old_quality;
// we have to roll back and reduce the step size
for (Uint local_nd_id = 0 ; local_nd_id < nb_free_nodes ; ++local_nd_id) {
Uint global_nd_id = this->free_nodes[local_nd_id];
for (Uint i = 0 ; i < DIM ; ++i) {
this->nodes[global_nd_id][i] += this->grad[DIM*local_nd_id+i];
}
}
initial_step_size *= .5;
keep_max(initial_step_size, min_step_size);
if (verbose) {
std::cout << "Reduced step size to " << initial_step_size << endl;
}
active_max_iter = max_iter + iter;
}
}
if (verbose) {
std::cout << " Success: reached a quality of " << quality << " in " << iter
<< " iteration";
if (iter != 1){
std::cout << "s";
}
std::cout << std::endl;
}
return iter;
}
/* -------------------------------------------------------------------------- */
template class MiniMesh<twoD>;
template class MiniMesh<threeD>;
c4science · Help