// -*- C++ -*-
// $RCSfile: fiddlemesh.C,v $
// $Revision: 1.15 $
// $Author: langer $
// $Date: 2004/10/23 00:48:47 $
/* This software was produced by NIST, an agency of the U.S. government,
* and by statute is not subject to copyright in the United States.
* Recipients of this software assume all responsibilities associated
* with its operation, modification and maintenance. However, to
* facilitate maintenance we ask that before distributing modifed
* versions of this software, you first contact the authors at
* oof_manager@ctcms.nist.gov.
*/
// Stuff for adjusting the node positions.
#include "adaptmesh.h"
#include "amtriangle.h"
#include "amtriangleiterator.h"
#include "fiddlemesh.h"
#include "goof.h"
#include "material.h"
#include "menuinterrupt.h"
#include "meshcmds.h"
#include "random.h"
#include "timer.h"
#include "stdlib.h"
#include <math.h>
static const double minarea = 0.01; // minimum allowed area of a triangle
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// find all the different types of pixels
void Goof::categorize_pixels() {
// Compare time stamps. Don't do anything unless groups and
// materials have changed since the last time the pixels were
// categorized.
if(pixels_categorized > groups_changed &&
pixels_categorized > materials_changed)
return;
++pixels_categorized;
weed_all_groups();
// keep a list of pixels with unique properties
catalog.resize(0);
// loop over pixels
Cell_coordinate pixel;
ArrayIterator iter(material);
while(iter(pixel)) {
// does this pixel fit an existing category?
int fits = 0;
for(int i=0; i<catalog.capacity() && !fits; i++) {
if(material[pixel] == material[catalog[i]] &&
pixelgrouplist[pixel] == pixelgrouplist[catalog[i]]) { // yes
pixelcategory[pixel] = i;
fits = 1;
}
}
if(!fits) { // a new category of pixel
pixelcategory[pixel] = catalog.capacity();
catalog.grow(1, pixel);
}
}
Ncategories = catalog.capacity();
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// The sum of E over all triangles is the function to be minimized by
// fiddling with the node positions. E has two parts: a homogeneity
// energy that prefers homogenous triangles, and an angle energy that
// prefers fat triangles with no sharp angles.
// The area of all pixels of type i within the triangle is A*a[i],
// where A is the area of the triangle. The homogeneiety energy
// function is
// E_h = \product_i (1-a[i])*normalization
// The energy is minimized (E_h=0) if the triangle is homogeneous, so
// that one a[i] is one and the rest are zero. It it normalized so
// that its maximum value (obtained when all a[i]=1/N) is one.
double AdaptiveMesh::alpha = 0.33;
double AMTriangle::E() {
if(AdaptiveMesh::alpha == 0.0)
return Eangle();
return AdaptiveMesh::alpha*Ehomogeneity() + (1-AdaptiveMesh::alpha)*Eangle();
}
double AMTriangle::Ehomogeneity() {
// is the saved E current?
if(mesh->goof->materials_changed < Ecalculated
&& node[0]->nodemoved < Ecalculated
&& node[1]->nodemoved < Ecalculated
&& node[2]->nodemoved < Ecalculated)
return currentE;
mesh->goof->categorize_pixels(); // find all pixel types
int N = mesh->goof->Ncategories;
if(N <= 1) return 0;
Vec<double> ac(N, 0.0); // area of each type
const Array<int> &pixelcategory = mesh->goof->pixelcategory;
for(AMTriangleIterator ti(*this); !ti.end(); ++ti) {
const Cell_coordinate pixel = (*this)[ti];
double ai = intersection(pixel);
ac[pixelcategory[pixel]] += ai;
}
double e = 1.0;
double a = area();
double invarea = 1./a;
double norm = N/(N-1.0); // 1/(1-1/N)
int j;
for(j=0; j<N; j++) {
ac[j] *= invarea; // convert to fractional areas
e *= (1-ac[j])*norm;
}
lastE = currentE;
currentE = e;
lastEcalculated = Ecalculated;
++Ecalculated;
return e;
}
void AMTriangle::revertE() {
Ecalculated = lastEcalculated;
currentE = lastE;
}
void AMNode::revert() {
move_back();
for(int i=0; i<triangle.capacity(); i++)
triangle[i]->revertE();
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// Energy function that is a minimum when the triangle has no small angles.
double AMTriangle::Eangle() const {
static const double norm = 36./sqrt(3.); // L^2/A for an equilateral triangle
double L = 0;
for(int i=0; i<3; i++) {
MeshCoord edge = node[i]->coord() - node[(i+1)%3]->coord();
L += edge.norm();
}
return 1.0 - area()/(L*L)*norm;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
double AMNode::E() const {
double e = 0;
for(int i=0; i<triangle.capacity(); i++)
e += triangle[i]->E();
return e;
}
double AMNode::Eangle() const {
double e = 0;
for(int i=0; i<triangle.capacity(); i++)
e += triangle[i]->Eangle();
return e;
}
double AMNode::Ehomogeneity() const {
double e = 0;
for(int i=0; i<triangle.capacity(); i++)
e += triangle[i]->Ehomogeneity();
return e;
}
double AdaptiveMesh::E() {
double e = 0;
for(AMIterator it(this, AMI_ALL); !it.end(); ++it)
e += (*this)[it]->E();
return e;
}
void AdaptiveMesh::Erange(double &min, double &max) {
min = 2;
max = -1;
for(AMIterator it(this, AMI_ALL); !it.end(); ++it) {
if((*this)[it]->active()) {
double e = (*this)[it]->E();
if(e > max) max = e;
if(e < min) min = e;
}
}
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
int AMNode::MCmove(double T, double delta, double *deltaE) {
// Choose a random displacement with a Gaussian distribution, but
// don't move node off of the boundaries.
double dx, dy;
if(lft() || rgt())
dx = 0;
else
dx = gasdev()*delta;
if(top() || btm())
dy = 0;
else
dy = gasdev()*delta;
MeshCoord disp(dx, dy);
double estart = E(); // original energy
move_by(disp); // move to new position
// Make sure that topology hasn't changed!
for(int j=0; j<triangle.capacity(); j++)
if(triangle[j]->area() < 0.0) {
revert(); // oops. Move back to original position
return 0; // reject move
}
++nodemoved;
double efinal = E();
*deltaE = efinal - estart;
if(*deltaE < 0 || (T > 0 && exp(-*deltaE/T) > rndm())) // accept the move
return 1;
else { // don't accept the move
revert(); // move node back to where it was
return 0;
}
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
int AMNode::Laplacemove(double *deltaE) {
// find average position of neighboring nodes
MeshCoord avg;
int nnbr = 0;
for(int i=0; i<triangle.capacity(); i++) {
for(int j=0; j<3; j++)
if(triangle[i]->node[j] != this) {
avg += triangle[i]->node[j]->coord();
nnbr++;
}
}
avg /= nnbr;
// don't move off of boundaries
if(lft() || rgt())
avg.x = coord().x;
if(top() || btm())
avg.y = coord().y;
double estart = E();
move_to(avg);
double efinal = E();
*deltaE = efinal - estart;
if(*deltaE <= 0)
return 1;
// reject moves that increase E
revert();
return 0;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
Vec<AMNode*> AdaptiveMesh::activenodes() const {
Vec<AMNode*> anodes;
anodes.setphysicalsize(nodes.capacity());
for(int i=0; i<nodes.capacity(); i++)
if(nodes[i]->active())
anodes.grow(1, nodes[i]);
return anodes;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// create a list of active nodes in random order
static Vec<AMNode*> shufflenodes(const Vec<AMNode*> &nodelist) {
Vec<AMNode*> randomnode(nodelist);
int N = randomnode.capacity();
for(int i=0; i<N; i++) {
int x = int(i + rndm()*(N-i)); // pick a remaining node at random
if(x != i) { // swap random node with node i
AMNode *temp = randomnode[i];
randomnode[i] = randomnode[x];
randomnode[x] = temp;
}
}
return randomnode;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
void AdaptiveMesh::fiddleMC(double T, double delta,
const Vec<AMNode*> &nodelist)
{
Vec<AMNode*> randomnode(shufflenodes(nodelist));
double deltaE = 0; // cumulative change in energy
// try to move each node
int accepted = 0;
int rejected = 0;
for(int i=0; i<randomnode.capacity(); i++) { // loop over shuffled nodes
double dE;
if(randomnode[i]->MCmove(T, delta, &dE)) {
deltaE += dE;
accepted++;
// current_goof->redraw();
}
else
rejected++;
}
materials_need_recomputing();
pixels_listed = 0; // force update of pixel lists in each triangle
garcon()->msout << ms_info << "dE = " << deltaE << " acceptance rate = "
<< float(accepted)/(accepted + rejected) << endl
<< ms_normal;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
void AdaptiveMesh::fiddleLaplace(const Vec<AMNode*> &nodelist) {
Vec<AMNode*> randomnode(shufflenodes(nodelist));
double deltaE = 0;
int accepted = 0;
int rejected = 0;
for(int i=0; i<randomnode.capacity(); i++) {
double dE;
if(randomnode[i]->Laplacemove(&dE)) {
deltaE += dE;
accepted++;
}
else
rejected++;
}
materials_need_recomputing();
pixels_listed = 0;
garcon()->msout << ms_info << "dE = " << deltaE << " acceptance rate = "
<< float(accepted)/(accepted + rejected)
<< endl << ms_normal;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
int AdaptiveMesh::fiddleDownhill(double delta) {
static const double epsilon = 1.e-5;
MeshCoord dx(epsilon, 0);
MeshCoord dy(0, epsilon);
Vec<MeshCoord> gradE(nodes.capacity());
Interrupter stop;
double e0 = E();
double gradmax = 0; // maximum gradient found
for(int i=0; i<nodes.capacity() && !stop(); i++) {
double e0 = nodes[i]->E();
double f=1;
// compute dE/dx
if(!(nodes[i]->rgt() || nodes[i]->lft())) {
nodes[i]->move_by(dx); // move the node
while(!nodes[i]->areas_ok(minarea)) { // check for negative area
nodes[i]->revert(); // move back
f *= -0.5; // try again with half the motion in
nodes[i]->move_by(f*dx); // the opposite direction
}
gradE[i].x = (nodes[i]->E() - e0)/(f*epsilon);
nodes[i]->revert(); // move back
}
// compute dE/dy
if(!(nodes[i]->top() || nodes[i]->btm())) {
f=1;
nodes[i]->move_by(dy); // move the node
while(!nodes[i]->areas_ok(minarea)) { // check for negative area
nodes[i]->revert(); // move back
f *= -0.5; // try again with half the motion in
nodes[i]->move_by(f*dy); // the opposite direction
}
gradE[i].y = (nodes[i]->E() - e0)/(f*epsilon);
nodes[i]->revert(); // move back
}
double g = dot(gradE[i], gradE[i]);
if(g > gradmax) gradmax = g;
}
// if(stop()) return;
gradmax = sqrt(gradmax);
// move all nodes
double move = delta/gradmax;
for(int j=0; j<nodes.capacity(); j++)
nodes[j]->move_by(-move*gradE[j]);
if(!areas_ok(minarea) || E() > e0) { // reject move!
// move all nodes back
for(int k=0; k<nodes.capacity(); k++)
nodes[k]->revert();
return 0;
}
garcon()->msout << ms_info << "dE = " << E() - e0
<< " max grad = " << gradmax << endl << ms_normal;
materials_need_recomputing();
pixels_listed = 0;
return 1;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// node3 2 1
// /\ /|\
// Change / \ to / | \
// / \ / | \
// / A \ / | \
// node0 /________\ node2 0/ | \0
// \ / \ A' | B' /
// \ B / \ | /
// \ / \ | /
// \ / \ | /
// \/node1 1\|/2 <--- nodenumbers of new triangles
//
int AdaptiveMesh::swap_edges(AMTriangle *A, AMTriangle *B,
AMTriangle *&Anew, AMTriangle *&Bnew)
{
int i;
// identify the nodes in common
int nbr_a = 0; // which neighbor of A is B?
int nbr_b = 0; // which neighbor of B is A?
while(nbr_a < 3 && A->neighbor[nbr_a] != B) ++nbr_a;
while(nbr_b < 3 && B->neighbor[nbr_b] != A) ++nbr_b;
if(nbr_a == 3 || nbr_b == 3) {
cerr << "AdaptiveMesh::swap_edges: bad pair!" << endl;
return 0;
}
AMNode *node0 = A->node[(nbr_a+1)%3];
AMNode *node1 = B->node[nbr_b];
AMNode *node2 = B->node[(nbr_b+1)%3];
AMNode *node3 = A->node[nbr_a];
// check that resulting triangles will have positive areas!
if((trianglearea(node0->coord(), node1->coord(), node3->coord()) <= minarea)||
(trianglearea(node2->coord(), node3->coord(), node1->coord()) <= minarea))
return 0;
// remove old triangles from lists in the nodes
for(i=0; i<3; i++) {
A->node[i]->remove_triangle(A);
B->node[i]->remove_triangle(B);
}
// create new triangles
Anew = new AMTriangle(A->parent, node0, node1, node3);
Bnew = new AMTriangle(B->parent, node2, node3, node1);
Anew->mesh = this;
Bnew->mesh = this;
// update neighbor lists
int nbr_a1 = (nbr_a+1)%3;
int nbr_a2 = (nbr_a+2)%3;
int nbr_b1 = (nbr_b+1)%3;
int nbr_b2 = (nbr_b+2)%3;
Anew->neighbor[0] = Bnew;
Anew->neighbor[1] = A->neighbor[nbr_a2];
Anew->neighbor[2] = B->neighbor[nbr_b1];
Bnew->neighbor[0] = Anew;
Bnew->neighbor[1] = B->neighbor[nbr_b2];
Bnew->neighbor[2] = A->neighbor[nbr_a1];
if(A->neighbor[nbr_a1])
A->neighbor[nbr_a1]->replace_neighbor(A, Bnew);
if(A->neighbor[nbr_a2])
A->neighbor[nbr_a2]->replace_neighbor(A, Anew);
if(B->neighbor[nbr_b1])
B->neighbor[nbr_b1]->replace_neighbor(B, Anew);
if(B->neighbor[nbr_b2])
B->neighbor[nbr_b2]->replace_neighbor(B, Bnew);
// update parents or root list
if(A->parent) {
if(A->parent->child[0] == A)
A->parent->child[0] = Anew;
else
A->parent->child[1] = Anew;
}
else { // A is a root triangle
for(int k=0; k<nroots; k++)
if(root[k] == A) {
root[k] = Anew;
break;
}
}
if(B->parent) {
if(B->parent->child[0] == B)
B->parent->child[0] = Bnew;
else
B->parent->child[1] = Bnew;
}
else { // B is a root triangle
for(int k=0; k<nroots; k++)
if(root[k] == B) {
root[k] = Bnew;
break;
}
}
// make sure all derivative quantities will be recalculated
Anew->representative_material_time.backdate();
Bnew->representative_material_time.backdate();
Anew->group_cell_time.backdate();
Bnew->group_cell_time.backdate();
Anew->meshgroups_time.backdate();
Bnew->meshgroups_time.backdate();
// Set the new triangles' inhibit_inheritance flags. If the
// triangles have the same material and both flags are set, then it
// makes sense to set the flags in the new triangles. Otherwise,
// both new triangles should have the flags unset.
if(A->material() && B->material() && *A->material() == *B->material() &&
A->inhibit_inheritance && B->inhibit_inheritance)
{
Anew->set_material(A->material(), true/* inhibit automatic assignment*/);
Bnew->set_material(B->material(), true/* inhibit automatic assignment*/);
}
else {
Anew->inhibit_inheritance = false;
Bnew->inhibit_inheritance = false;
}
return 1;
}
void AdaptiveMesh::unswap_edges(AMTriangle *old1, AMTriangle *old2,
AMTriangle *new1, AMTriangle *new2)
{
// Replace new1 and new2 with old1 and old2. All links in the old
// triangles are still ok, but the links from other triangles have
// to be restored.
int i;
for(i=0; i<3; i++) {
// restore node indices
new1->node[i]->remove_triangle(new1);
new2->node[i]->remove_triangle(new2);
old1->node[i]->add_triangle(old1);
old2->node[i]->add_triangle(old2);
}
if(!new1->parent) root.remove(new1);
if(!new2->parent) root.remove(new2);
// restore links in parents
if(old1->parent) {
if(old1->parent->child[0] == new1 || old1->parent->child[0] == new2)
old1->parent->child[0] = old1;
else
old1->parent->child[1] = old1;
}
else // old1 is a root triangle
root.grow(1, old1);
if(old2->parent) {
if(old2->parent->child[0] == new1 || old2->parent->child[0] == new2)
old2->parent->child[0] = old2;
else
old2->parent->child[1] = old2;
}
else // old2 is a root triangle
root.grow(1, old2);
// restore links in neighbors
for(i=0; i<3; i++) {
if(old1->neighbor[i] && old1->neighbor[i] != old2)
if(!old1->neighbor[i]->replace_neighbor(new1, old1))
old1->neighbor[i]->replace_neighbor(new2, old1);
if(old2->neighbor[i] && old2->neighbor[i] != old1)
if(!old2->neighbor[i]->replace_neighbor(new1, old2))
old2->neighbor[i]->replace_neighbor(new2, old2);
}
// force derivative quantities to be updated
old1->representative_material_time.backdate();
old2->representative_material_time.backdate();
old1->group_cell_time.backdate();
old2->group_cell_time.backdate();
old1->meshgroups_time.backdate();
old2->meshgroups_time.backdate();
}
void AdaptiveMesh::test_swap() {
// find selected triangles
AMTriangle *t1 = 0;
AMTriangle *t2 = 0;
for(AMIterator i(this, AMI_SELECTED); !i.end(); ++i)
if(!t1)
t1 = (*this)[i];
else if(!t2)
t2 = (*this)[i];
else {
garcon()->msout << ms_error << "Please select only two triangles!"
<< endl << ms_normal;
return;
}
if(!t1 || !t2) {
garcon()->msout << ms_error << "Please select two triangles!"
<< endl << ms_normal;
return;
}
double e0 = t1->E() + t2->E();
AMTriangle *new1, *new2;
if(swap_edges(t1, t2, new1, new2)) {
goof->triangle_destroyed(t1);
goof->triangle_destroyed(t2);
delete t1;
delete t2;
new1->select();
new2->select();
materials_need_recomputing();
inherit_pixel_materials();
groups_need_recomputing();
inherit_pixel_groups();
double e1 = new1->E() + new2->E();
garcon()->msout << ms_info << "swap succeeded! dE=" << (e1 - e0)
<< endl << ms_normal;
current_goof->redraw();
}
else
garcon()->msout << ms_info << "swap failed!" << endl << ms_normal;
}
static int compareT(const void *a, const void *b) {
double ea = (*(AMTriangle**) a)->E();
double eb = (*(AMTriangle**) b)->E();
if(ea > eb) return -1;
if(ea < eb) return 1;
return 0;
}
void AdaptiveMesh::swap_worst() {
// make a list of all triangles, sorted by decreasing energy
Vec<AMTriangle*> tangleangle; // Naomi's word for "triangle"
tangleangle.setphysicalsize(ntriangles());
for(AMIterator it(this, AMI_ALL); !it.end(); ++it) {
if((*this)[it]->active()) {
tangleangle.grow(1, (*this)[it]); // unsorted list of active triangles
(*this)[it]->swapped = 0; // reset flag
}
}
::qsort(&tangleangle[0], tangleangle.capacity(), sizeof(AMTriangle*),
compareT); // sort it
// Triangles that have been replaced can't be deleted immediately,
// since the tangleangle list still points to them. So keep a list
// of defunct triangles and delete them at the end.
Vec<AMTriangle*> defunct(0, BlockSize(10));
// look for pairs to swap
int nswapped = 0;
double dE = 0;
for(int i=0; i<tangleangle.capacity(); i++) {
if(!tangleangle[i]->swapped) {
AMTriangle *tri0 = tangleangle[i]; // current triangle is tri0
// Mark this triangle as swapped, even if nothing is done. If no
// swap is performed now, it won't be done when a neighbor is
// examined either, so marking it now saves time later.
tri0->swapped = 1;
// find neighbor (tri1) with highest energy
AMTriangle *tri1 = 0;
double en = -1.e30; // smaller than any possible E
for(int k=0; k<3; k++) {
AMTriangle *nbr = tri0->neighbor[k];
if(nbr && !nbr->swapped) {
if(nbr->E() > en) {
tri1 = nbr;
en = nbr->E();
}
}
}
if(tri1) {
// try to swap tri0 and tri1. keep the swap if the energy decreases
double e0 = tri0->E() + tri1->E();
AMTriangle *new0, *new1;
if(swap_edges(tri0, tri1, new0, new1)) {
double e1 = new0->E() + new1->E();
if(e1 < e0) { // swap succeeded in lowering energy
tri0->swapped = 1;
tri1->swapped = 1;
new0->swapped = 1;
new1->swapped = 1;
defunct.grow(1, tri0);
defunct.grow(1, tri1);
nswapped++;
dE += e1 - e0;
}
else { // swap failed to lower the energy
unswap_edges(tri0, tri1, new0, new1);
delete new0;
delete new1;
}
}
}
}
}
// finally delete old (swapped) triangles
for(int j=0; j<defunct.capacity(); j++) {
goof->triangle_destroyed(defunct[j]);
delete defunct[j];
}
garcon()->msout << ms_info << "nswapped=" << nswapped
<< " dE=" << dE << endl << ms_normal;
materials_need_recomputing();
inherit_pixel_materials();
groups_need_recomputing();
inherit_pixel_groups();
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
int AdaptiveMesh::areas_ok(double min) const {
for(AMIterator it(this, AMI_ALL); !it.end(); ++it)
if((*this)[it]->area() <= min)
return 0;
return 1;
}
int AMNode::areas_ok(double min) const {
for(int i=0; i<triangle.capacity(); i++)
if(triangle[i]->area() <= min)
return 0;
return 1;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
FiddleCmd::FiddleCmd()
: CommandM("anneal", "optimize mesh by simulated annealing"),
T(0),
delta(1.0),
iterations(1)
{
AddArgument(this, "T", T);
AddArgument(this, "delta", delta);
AddArgument(this, "iterations", iterations);
}
CommandFn FiddleCmd::func() {
if(meshexists()) {
Timer timer;
Interrupter stop;
current_goof->dup_mesh();
Vec<AMNode*> nodelist(current_goof->mesh()->activenodes());
for(int i=0; i<iterations && !stop(); i++) {
current_goof->mesh()->fiddleMC(T, delta, nodelist);
if(AdaptiveMesh::continuous_redraw)
current_goof->redraw();
}
if(!AdaptiveMesh::continuous_redraw)
current_goof->redraw();
garcon()->msout << ms_info << timer << endl << ms_normal;
}
}
CommandM *FiddleCmd::clone() const {
FiddleCmd *fc = new FiddleCmd;
fc->T = T;
fc->delta = delta;
fc->iterations = iterations;
return fc;
}
FiddleDownHillCmd::FiddleDownHillCmd()
: CommandM("relax", "optimize mesh by relaxation"),
delta(1),
iterations(1)
{
AddArgument(this, "delta", delta);
AddArgument(this, "iterations", iterations);
}
CommandFn FiddleDownHillCmd::func() {
if(meshexists()) {
Interrupter stop;
current_goof->dup_mesh();
for(int i=0; i<iterations && !stop(); i++) {
if(current_goof->mesh()->fiddleDownhill(delta)) {
if(AdaptiveMesh::continuous_redraw)
current_goof->redraw();
}
else {
garcon()->msout << ms_error
<< "Failed to move nodes. Try reducing delta."
<< endl << ms_normal;
break;
}
}
if(!AdaptiveMesh::continuous_redraw)
current_goof->redraw();
}
}
CommandM *FiddleDownHillCmd::clone() const {
FiddleDownHillCmd *fdhc = new FiddleDownHillCmd;
fdhc->delta = delta;
fdhc->iterations = iterations;
return fdhc;
}
FiddleLaplaceCmd::FiddleLaplaceCmd()
: CommandM("smooth", "optimize mesh by Laplacian smoothing"),
iterations(3)
{
AddArgument(this, "iterations", iterations);
}
CommandFn FiddleLaplaceCmd::func() {
if(meshexists()) {
Timer timer;
Interrupter stop;
current_goof->dup_mesh();
Vec<AMNode*> nodelist(current_goof->mesh()->activenodes());
for(int i=0; i<iterations && !stop(); i++) {
current_goof->mesh()->fiddleLaplace(nodelist);
if(AdaptiveMesh::continuous_redraw)
current_goof->redraw();
}
if(!AdaptiveMesh::continuous_redraw)
current_goof->redraw();
garcon()->msout << ms_info << timer << endl << ms_normal;
}
}
CommandM *FiddleLaplaceCmd::clone() const {
FiddleLaplaceCmd *flc = new FiddleLaplaceCmd;
flc->iterations = iterations;
return flc;
}
CommandFn swapedges() {
if(meshexists()) {
current_goof->dup_mesh();
current_goof->mesh()->test_swap();
current_goof->redraw();
}
}
CommandFn swapworst() {
if(meshexists()) {
current_goof->dup_mesh();
current_goof->mesh()->swap_worst();
current_goof->redraw();
}
}
syntax highlighted by Code2HTML, v. 0.9.1