// -*- C++ -*-
// $RCSfile: amtriangle.C,v $
// $Revision: 1.44 $
// $Author: langer $
// $Date: 2004/10/22 13:51:08 $
/* 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.
*/
//#include "config.h"
#include "adaptmesh.h"
#include "amtriangle.h"
#include "amtriangleiterator.h"
#include "colorutils.h"
#include "elector.h"
#include "goof.h"
#include "imagecanvas.h"
#include "material.h"
#include "meshcmds.h"
#include "pixelgroups.h"
#include "saveconfig.h"
#include "word.h"
#include <iostream.h>
#include <math.h>
#include "stdlib.h"
#ifdef HAVE_MALLOC_H
#include <malloc.h>
#endif // HAVE_MALLOC_H
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
const unsigned char AMNode::top_ = 1;
const unsigned char AMNode::btm_ = 2;
const unsigned char AMNode::lft_ = 4;
const unsigned char AMNode::rgt_ = 8;
// put a node between two nodes
AMNode::AMNode(AdaptiveMesh *mesh, const AMNode *n1, const AMNode *n2)
: coord_(0.5*(n1->coord() + n2->coord())),
edgeflags(n1->edgeflags & n2->edgeflags)
{
index = mesh->nodes.capacity();
mesh->nodes.grow(1, this);
}
AMNode::AMNode(AdaptiveMesh *mesh, double x, double y,
unsigned char e)
: coord_(x, y),
edgeflags(e)
{
index = mesh->nodes.capacity();
mesh->nodes.grow(1, this);
}
void AMNode::move_to(const MeshCoord &p) {
lastmoved = nodemoved;
lastcoord = coord_;
++nodemoved;
coord_ = p;
for(int i=0; i<triangle.capacity(); i++)
++triangle[i]->lastmod_time;
}
void AMNode::move_by(const MeshCoord &dp) {
lastmoved = nodemoved;
lastcoord = coord_;
++nodemoved;
coord_ += dp;
for(int i=0; i<triangle.capacity(); i++)
++triangle[i]->lastmod_time;
}
void AMNode::move_back() {
coord_ = lastcoord;
for(int i=0; i<triangle.capacity(); i++)
++triangle[i]->lastmod_time;
nodemoved = lastmoved;
}
AMNode *AMNode::copy() const {
AMNode *newnode = new AMNode;
newnode->index = index;
newnode->edgeflags = edgeflags;
newnode->coord_ = coord_;
newnode->lastcoord = lastcoord;
newnode->lastmoved = lastmoved;
// don't copy triangle list -- will be done when triangles are copied
return newnode;
}
bool AMNode::active() const {
return current_goof->active(coord_);
}
ostream &operator<<(ostream &os, const AMNode &node) {
return os << node.coord_;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// Routines for saving and reading nodes in configuration files
// Don't save the index, because nodes will be stored in order, and
// don't save the timestamps, so that the quantities they govern will
// be recomputed automatically.
void AMNode::save(ostream &file) const {
int ef = edgeflags; // writing unsigned chars is dangerous
file << ef << " " << coord_ << " " << lastcoord << endl;
}
void AMNode::read(AdaptiveMesh *mesh, istream &file) { // static
int edgeflags;
MeshCoord coord, lastcoord;
file >> edgeflags >> coord >> lastcoord;
// constructor puts node in list in mesh
if(file) {
AMNode *knowed = new AMNode(mesh, coord.x, coord.y, edgeflags);
knowed->lastcoord = lastcoord;
}
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// Comparison operators on Nodes are used in
// AMTriangle::list_other_nodes and when comparing AMTriangles. We'll
// never be comparing Nodes that are in different Meshes, so it's
// sufficient to compare their indices.
int operator!=(const AMNode &n1, const AMNode &n2) {
return (n1.index != n2.index);
}
int operator<(const AMNode &n1, const AMNode &n2) {
return (n1.index < n2.index);
}
int operator>(const AMNode &n1, const AMNode &n2) {
return (n1.index > n2.index);
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
Enum<TriangleDivisionMethod> AMTriangle::divisionmethod(TDM_LONGESTEDGE);
void AMTriangle::init() {
// call this just once
static int once = 0;
if(once) return;
Enum<TriangleDivisionMethod>::name(TDM_NEWESTNODE, "newest_node");
Enum<TriangleDivisionMethod>::name(TDM_LONGESTEDGE, "longest_edge");
Enum<TriangleDivisionMethod>::name(TDM_E, "smallest_E");
once = 1;
return;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
const unsigned char AMTriangle::SELECTED = 0x1;
// There used to be more flags defined here.
static int trianglecount = 0;
AMTriangle::AMTriangle(AMTriangle *mom,
AMNode *n0, AMNode *n1, AMNode *n2)
: parent(mom),
status_(0),
representative_material(0),
swapped(0),
inhibit_inheritance(0),
inhibit_groupinheritance(0),
index(trianglecount++) // used for debugging and reused for output
{
node[0] = n0;
node[1] = n1;
node[2] = n2;
n0->add_triangle(this);
n1->add_triangle(this);
n2->add_triangle(this);
child[0] = child[1] = 0;
if(parent) {
mesh = parent->mesh;
if(parent->selected()) select();
generation = parent->generation+1;
initial_generation = parent->initial_generation;
if(mesh->depth < generation) mesh->depth = generation;
}
else {
mesh = 0; // hopefully will be set later
generation = 0;
initial_generation = 0;
}
neighbor[0] = neighbor[1] = neighbor[2] = 0;
Ecalculated.backdate(); // force initial calculation of E
representative_material_time.backdate();
++lastmod_time;
}
AMTriangle::~AMTriangle() {
// Don't call AMNode::remove_triangle(this), because when the mesh
// is destructed, the nodes are destructed before the triangles!
mesh->goof->triangle_destroyed(this);
if(child[0]) delete child[0];
if(child[1]) delete child[1];
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
void AMTriangle::save(ostream &file) const {
file << node[0]->index << " "
<< node[1]->index << " "
<< node[2]->index << " "
<< selected() << " "
<< (child[0] != 0) << endl;
// added for config_version 3 -- all triangles store their own materials
file << inhibit_inheritance
<< " " << representative_material->index() << endl;
// added for config_version 3 -- all triangles store their own groups
file << inhibit_groupinheritance
<< " " << meshgroups.capacity();
for(int i=0; i<meshgroups.capacity(); i++) {
file << " " << meshgroups[i]->query_name();
}
file << endl;
if(child[0]) {
child[0]->save(file);
child[1]->save(file);
}
}
AMTriangle *AMTriangle::read(istream &file, AdaptiveMesh *mesh,
AMTriangle *parent)
{
int nodeindex[3];
file >> nodeindex[0] >> nodeindex[1] >> nodeindex[2];
if(!file) return 0;
bool slctd;
file >> slctd;
if(!file) return 0;
bool is_parent;
file >> is_parent;
if(!file) return 0;
bool inhibited = 0;
int mat = 0;
if(configversion() == 2) {
file >> inhibited;
if(!file) return 0;
if(inhibited) {
file >> mat;
if(!file) return 0;
}
}
bool group_inhibited = 0;
Vec<Word*> glist;
if(configversion() >= 3) {
file >> inhibited;
if(!file) return 0;
file >> mat;
if(!file) return 0;
file >> group_inhibited;
if(!file) return 0;
int ngrps;
file >> ngrps;
if(!file) return 0;
for(int i=0; i<ngrps; i++) {
Word *grp = new Word;
file >> *grp;
if(!file) return 0;
glist.grow(1, grp);
}
}
AMTriangle *triangle = new AMTriangle(parent,
mesh->nodes[nodeindex[0]],
mesh->nodes[nodeindex[1]],
mesh->nodes[nodeindex[2]]);
if(slctd)
triangle->select();
triangle->mesh = mesh;
triangle->set_material(Material::allmaterials[mat], inhibited);
if(group_inhibited)
triangle->inhibit_group_inheritance();
for(int i=0; i<glist.capacity(); i++) {
mesh->find_group(*glist[i])->append(triangle);
delete glist[i];
}
if(is_parent) {
// Oops. AMTriangle constructor called Node::add_triangle(), but
// we're reading a refined triangle. When triangles are refined,
// AMTriangle::divide() calls Node::remove_triangle(). We need to
// call that by hand here. This is ugly, and should be done
// better in some later version.
mesh->nodes[nodeindex[0]]->remove_triangle(triangle);
mesh->nodes[nodeindex[1]]->remove_triangle(triangle);
mesh->nodes[nodeindex[2]]->remove_triangle(triangle);
// read the children
triangle->child[0] = read(file, mesh, triangle);
triangle->child[1] = read(file, mesh, triangle);
}
return triangle;
}
ostream &operator<<(ostream &os, const AMTriangle &triangle) {
os << "triangle "
<< *triangle.node[0] << " "
<< *triangle.node[1] << " "
<< *triangle.node[2];
return os;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
void AMTriangle::mark_for_division() {
// See the comment in AMTriangle::divide() wrt the order in which
// triangles are marked for division.
mesh->dividelist.grow(1, this);
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
void AMTriangle::divide() {
if(divided())
return;
// . <-- node div of parent
// /|\
// / \
// *this ----> / | \
// / 0 1 \ <--- child numbers
// / | \
// node1 of parent -> /_____._____\ <--- node2 of parent
// \ ^newnode
// \ /
// \ nbr /
// \ /
//
// .
// /|\
// Node numbering /1|2\ For all triangles, neighbor[i]
// of children / | \ is opposite node[i]
// / | \
// / | \
// /2___0|0___1\
bool divide_child[2];
bool divide_nbr = false;
divide_child[0] = divide_child[1] = false;
int div; // node through which to divide triangle
if((mesh->max_divisions > 0 && division_depth() >= mesh->max_divisions)
|| (area() <= mesh->min_area)) {
// Triangle has been divided too much already. Divide again only if
// division is forced by a neighbor, and divide compatibly with
// the neighbor, so that no further divisions will be needed. If
// two or more neighbors force division, then this should use the
// current divisionmethod to chose between them, but it doesn't.
bool forced = false;
for(int i=0; i<3 && !forced; i++)
if(neighbor[i] && neighbor[i]->divided()) {
// If the neighbor has been divided, then it has to have been
// divided through the edge it shares with this triangle. If
// it had been divided in some other direction, then one of
// its children would be this triangle's neighbor instead!
div = i;
forced = true;
}
if(!forced) // division wasn't forced
return; // don't divide
}
else { // triangle hasn't been divided too much
div = divider(); // choose optimal direction
}
divided_node = div;
AMTriangle *nbr = neighbor[div]; // neighbor across divided edge
AMNode *newnode = 0; // new node
if(nbr && nbr->divided())
newnode = nbr->child[0]->node[0];
else
newnode = new AMNode(mesh, node[(div+1)%3], node[(div+2)%3]);
// create children
int node1 = (div + 1) % 3;
int node2 = (div + 2) % 3;
child[0] = new AMTriangle(this, newnode, node[div], node[node1]);
child[1] = new AMTriangle(this, newnode, node[node2], node[div]);
// update each node's list of triangles
node[0]->remove_triangle(this);
node[1]->remove_triangle(this);
node[2]->remove_triangle(this);
// update the neighbor pointers
child[0]->neighbor[2] = child[1];// children are neighbors of each other
child[1]->neighbor[1] = child[0];
child[0]->neighbor[0] = neighbor[node2]; // neighbors of parent are
child[1]->neighbor[0] = neighbor[node1]; // neighbors of the children
if(neighbor[node1])
neighbor[node1]->replace_neighbor(this, child[1]);
if(neighbor[node2])
neighbor[node2]->replace_neighbor(this, child[0]);
// Mark children for division if they have neighbors that have been
// divided.
if(neighbor[node2] && neighbor[node2]->divided())
divide_child[0] = true;
if(neighbor[node1] && neighbor[node1]->divided())
divide_child[1] = true;
if(nbr) {
if(nbr->divided()) {
// cerr << "Neighbor " << nbr->index << " has been divided" << endl;
// Neighbor has been divided compatibly, so its descendents are
// the neighbors of this triangle's children. If it's been
// divided incompatibly, then it shouldn't be a neighbor-- one
// of its descendents should be the neighbor!
// The neighbors children could have already been divided, so
// they may not be the neighbors of this triangle's
// children. Find which of the descendents of nbr are adjacent
// to this triangle's children.
AMTriangle *nbr0 =
nbr->child[1]->find_nbr_child(this, newnode, node[node1]);
child[0]->neighbor[1] = nbr0;
if(nbr0) {
nbr0->neighbor[nbr0->neighbor_no(newnode, node[node1])] = child[0];
if(nbr0->divided())
// nbr0 is divided through the edge it shares with child[0]
child[0]->mark_for_division(); // should be divide_child[0] = true; ?
}
AMTriangle *nbr1 =
nbr->child[0]->find_nbr_child(this, newnode, node[node2]);
child[1]->neighbor[2] = nbr1;
if(nbr1) {
nbr1->neighbor[nbr1->neighbor_no(newnode, node[node2])] = child[1];
if(nbr1->divided())
child[1]->mark_for_division();
}
}
else {
divide_nbr = true;
child[0]->neighbor[1] = 0; // will be assigned when nbr is divided
child[1]->neighbor[2] = 0;
}
}
for(int i=0; i<2; i++)
if(divide_child[i])
child[i]->mark_for_division();
if(divide_nbr)
nbr->mark_for_division();
// transfer groups to children
child[0]->meshgroups.resize(meshgroups.capacity());
child[1]->meshgroups.resize(meshgroups.capacity());
child[0]->inhibit_groupinheritance = inhibit_groupinheritance;
child[1]->inhibit_groupinheritance = inhibit_groupinheritance;
for(int i=0; i<meshgroups.capacity(); i++) {
child[0]->meshgroups[i] = meshgroups[i];
child[1]->meshgroups[i] = meshgroups[i];
meshgroups[i]->replace(this, child[0]);
meshgroups[i]->append(child[1]);
}
meshgroups.smash();
// transfer material to children
child[0]->representative_material = representative_material;
child[1]->representative_material = representative_material;
child[0]->inhibit_inheritance = inhibit_inheritance;
child[1]->inhibit_inheritance = inhibit_inheritance;
// make sure material is updated as soon as possible:
child[0]->representative_material_time.backdate();
child[1]->representative_material_time.backdate();
// parent is no longer selected
unselect();
mesh->goof->triangle_destroyed(this);
#ifdef DEBUG
if(AdaptiveMesh::continuous_redraw)
mesh->goof->redraw();
#endif
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
void AMTriangle::listnbrs(ostream &os) const { // used for debugging
for(int i=0; i<3; i++)
if(neighbor[i])
os << neighbor[i]->index << " ";
else
os << "none ";
}
bool AMTriangle::replace_neighbor(AMTriangle *oldnbr, AMTriangle *newnbr) {
for(int i=0; i<3; i++)
if(neighbor[i] == oldnbr) {
neighbor[i] = newnbr;
return true; // replacement was successful
}
return false; // replacement wasn't made
}
// Find which neighbor contains nodes n1 and n2, assuming that this
// triangle actually has nodes n1 and n2.
int AMTriangle::neighbor_no(const AMNode *n1, const AMNode *n2) const {
for(int i=0; i<3; i++) {
if(node[i] != n1 && node[i] != n2)
return i;
}
cerr << "Error in AMTriangle::neighbor_no()!" << endl;
return -1; // should throw an exception?
}
// Find the smallest child of this triangle that contains nodes n1 and
// n2. It is assumed that this and nbr both share nodes n1 and n2.
AMTriangle *AMTriangle::find_nbr_child(AMTriangle *nbr,
const AMNode *n1, const AMNode *n2)
{
// check that this contains n1 and n2
bool n1found = false;
bool n2found = false;
for(int i=0; i<3; i++) {
if(node[i] == n1) n1found = true;
if(node[i] == n2) n2found = true;
}
if(!n1found || !n2found) // got the wrong triangle
return 0;
if(!divided()) { // this is the neighbor
return this;
}
else {
AMTriangle *c0 = child[0]->find_nbr_child(nbr, n1, n2);
if(c0) return c0; // child[0] is the smallest neighbor
AMTriangle *c1 = child[1]->find_nbr_child(nbr, n1, n2);
if(c1) return c1; // child[1] is the smallest neighbor
// neither child contains nodes n1 and n2; this is the neighbor
return this;
}
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// return the node which, when the triangle is divided through it,
// produces the smallest E
Vec<int> AMTriangle::smallest_E_nodes() const {
Vec<int> smallest;
smallest.setphysicalsize(3);
double e[3];
double min = e[0] = E_tentative_divide(0);
int which = 0;
for(int i=1; i<3; i++) {
e[i] = E_tentative_divide(i);
if(e[i] < min) {
min = e[i];
which = i;
}
}
smallest.grow(1, which);
for(int i=1; i<3; i++) {
int j = (which+i)%3;
if(e[j] == min)
smallest.grow(1, j);
}
return smallest;
}
// find E of the children if the triangle were divided through node nodeno
double AMTriangle::E_tentative_divide(int nodeno) const {
AMTriangle *child0 = new AMTriangle;
AMTriangle *child1 = new AMTriangle;
AMNode *newnode = new AMNode;
newnode->move_to(0.5*(node[(nodeno+1)%3]->coord() +
node[(nodeno+2)%3]->coord()));
child0->node[0] = node[nodeno];
child0->node[1] = node[(nodeno+1)%3];
child0->node[2] = newnode;
child1->node[0] = node[nodeno];
child1->node[1] = newnode;
child1->node[2] = node[(nodeno+2)%3];
child0->mesh = mesh;
child1->mesh = mesh;
child0->child[0] = 0;
child0->child[1] = 0;
child1->child[0] = 0;
child1->child[1] = 0;
double e = child0->E() + child1->E();
delete child0;
delete child1;
delete newnode;
return e;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
int AMTriangle::divider() const { // which edge should be divided?
if(divisionmethod == TDM_NEWESTNODE)
return 0; // newest node is always node 0
// Some methods can be ambiguous, and return a vector of possible directions.
Vec<int> choices;
if(divisionmethod == TDM_LONGESTEDGE)
choices = longestedges();
else if(divisionmethod == TDM_E)
choices = smallest_E_nodes();
if(choices.capacity() == 1)
return choices[0];
// Have to decide between more than one option.
// Chose to divide compatibly with a neighbor, if possible.
for(int i=0; i<choices.capacity(); i++) {
if(neighbor[i] && neighbor[i]->divided() &&
neighbor[i]->neighbor[neighbor[i]->divided_node] == this)
return i;
}
// No neighbors have been divided. At this point we could try to
// divide through an edge with no neighbors, but that's hard to
// detect, since neighbor[i]==0 if there's no neighbor, OR if this
// triangle has been more refined than its neighbor. So give up,
// and return the first choice.
return choices[0];
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
const Vec<MeshGroup*> &AMTriangle::mesh_groups() {
mesh->inherit_pixel_groups();
return meshgroups;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// A triangle is active if any of its nodes is active.
bool AMTriangle::active() const {
for(int i=0; i<3; i++)
if(node[i]->active()) return 1;
return 0;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
int AMTriangle::nnbrs_selected() const {
int n = 0;
for(int i=0; i<3; i++)
if(neighbor[i] && neighbor[i]->selected())
n++;
return n;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// find closest pixel from a list of pixels
const Cell_coordinate
&AMTriangle::closest_pixel(const Vec<Cell_coordinate>& plist) const
{
MeshCoord middle = center();
int closest = 0;
double d = sq_distance(middle, mesh->cellcenter(plist[0]));
for(int i=1; i<plist.capacity(); i++) {
double dd = sq_distance(middle, mesh->cellcenter(plist[i]));
if(dd < d) {
dd = d;
closest = i;
}
}
return plist[closest];
}
// find closest pixel in whole image (ie pixel containing center of triangle)
Cell_coordinate AMTriangle::closest_pixel() const {
MeshCoord middle = center();
return Cell_coordinate(middle.x, middle.y);
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
int AMTriangle::contains(const MeshCoord &pt) const {
// This can fail due to round-off error!
MeshCoord s0 = node[0]->coord() - pt;
MeshCoord s1 = node[1]->coord() - pt;
MeshCoord s2 = node[2]->coord() - pt;
return ((s0%s1 >= 0) && (s1%s2 >= 0) && (s2%s0 >= 0));
}
// int AMTriangle::contains(const Cell_coordinate &pt) const {
// return contains(mesh->cellcenter(pt));
// }
AMTriangle *AMTriangle::child_containing(const MeshCoord &pt)
{
if(!child[0]) return this; // fn can't be const, if it returns "this".
const MeshCoord &nexus = child[0]->node[0]->coord();
const MeshCoord septum = child[0]->node[1]->coord() - nexus;
// Using child[0]->node[1] instead of simply node[0] allows this to
// work even if this triangle was divided through a node other than
// node[0].
if(septum % (pt - nexus) > 0)
return child[0]->child_containing(pt);
else
return child[1]->child_containing(pt);
}
AMTriangle *AMTriangle::child_containing(const Cell_coordinate &pt) {
return child_containing(mesh->cellcenter(pt));
}
// Quick check to see if a triangle is outside a rectangle. A
// negative response does NOT mean that the triangle is inside, just
// that the quick check failed to prove that it is outside.
bool AMTriangle::outside(const Rectangle &rect) const {
double x[3], y[3];
for(int i=0; i<3; i++) {
x[i] = node[i]->coord().x;
y[i] = node[i]->coord().y;
}
bool left = true;
bool right = true;
bool above = true;
bool below = true;
for(int i=0; i<3; i++) {
left = left && x[i] < rect.ll.x;
right = right && x[i] > rect.ur.x;
below = below && y[i] < rect.ll.y;
above = above && y[i] > rect.ur.y;
}
return left || right || above || below;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// Find the neighbor in the direction of a given pt, by finding which
// edge of the triangle intersects the line between the center of the
// triangle and the pt.
AMTriangle *AMTriangle::neighbor_towards(const MeshCoord &pt) const {
int i;
MeshCoord cntr(center());
MeshCoord dest(pt - cntr); // vector from center to pt (destination)
MeshCoord r[3]; // vectors from center to node
// cross products of dest with vectors from center to nodes of triangle
double cp[3];
double dp[3]; // dot products
for(i=0; i<3; i++) {
r[i] = node[i]->coord() - cntr;
cp[i] = cross(dest, r[i]);
dp[i] = dot(dest, r[i]);
}
// Two of the cross products have the same sign, so their nodes are
// on the same side of the line between pt and center. The third
// cross product corresponds to a node that is at one end of the
// edge we're seeking. If one cross product is zero, then the other
// two must have opposite signs, and it doesn't matter which we
// pick. If two or more cross products are zero, something is
// wrong!
// Once one end is found, the other end is the node with the larger
// *normalized* dot product with dest. The neighbor we want to return
// is the neighbor whose number is the same as the node that ISN'T
// on the intersecting edge, because the neighbor number is the
// number of the node opposite it in the triangle.
for(i=0; i<3; i++) {
int j = (i == 2 ? 0 : i+1);
int k = (i == 0 ? 2 : i-1);
if(cp[j]*cp[k] > 0) { // cross products j and k have the same sign
dp[j] /= r[j].norm();
dp[k] /= r[k].norm();
if(dp[j] > dp[k])
return neighbor[k];
else if(dp[k] > dp[j])
return neighbor[j];
else {
// dp[j] == dp[k]. Corners j and k are on a line
// perpendicular to dest, and they're on the same side of the
// line. Which corner we want depends on whether they're
// upstream or downstream from i. Since the vectors all
// originate at the center of the triangle, if dp[i]==0 then
// dp[j]!=dp[k], and so we don't have to worry about dp[i]==0
// here.
if(dp[i] < 0) {
// dp[j]==dp[k] > 0, and the point with the smaller cross
// product is the other end of the intersecting edge.
if(cp[j] < cp[k])
return neighbor[k];
else
return neighbor[j];
}
else { // dp[i] > 0
if(cp[j] > cp[k])
return neighbor[k];
else
return neighbor[j];
}
}
}
}
// If we got this far, then no pair of cross products have the same
// sign, so one of the corners must lie on the line between center
// and pt.
for(i=0; i<3; i++) {
if(cp[i] == 0.0) {
if(dp[i] > 0)
return neighbor[i==2? 0 : i+1]; // node i is upstream, choose other edge
else
return neighbor[i]; // node i is downstream
}
}
return 0;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
static const double third = 1./3.;
MeshCoord AMTriangle::center() const {
return third*(node[0]->coord() + node[1]->coord() + node[2]->coord());
}
double AMTriangle::area() const {
return trianglearea(node[0]->coord(), node[1]->coord(), node[2]->coord());
}
Vec<int> AMTriangle::longestedges() const {
Vec<int> longest;
longest.setphysicalsize(3);
double length[3];
double max = -1;
int which;
for(int i=0; i<3; i++) {
MeshCoord edge = node[(i+1)%3]->coord() - node[(i+2)%3]->coord();
length[i] = dot(edge, edge);
if(length[i] > max) {
max = length[i];
which = i;
}
}
longest.grow(1, which);
for(int i=1; i<3; i++) {
int j = (which+1)%3;
if(length[j] == max)
longest.grow(1, j);
}
// if(longest.capacity() > 1) {
// cerr <<"lengths = ";
// for(int i=0; i<3; i++)
// cerr << sqrt(length[i]) << " ";
// cerr << "longestedges found " << longest.capacity()
// << " edges with length " << sqrt(max) << endl;
// }
return longest;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
void AMTriangle::list_other_nodes(const AMNode *n, Vec<AMNode*> &other) const {
for(int i=0; i<3; i++)
if(node[i] != n) {
int found = 0;
for(int k=0; k<other.capacity() && !found; k++)
if(other[k] == node[i])
found = 1;
if(!found)
other.grow(1, node[i]);
}
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
void AMTriangle::draw(ImageCanvas &canvas) const {
XPoint xp[4];
xp[0] = canvas.mesh2xpoint(node[0]->coord());
xp[1] = canvas.mesh2xpoint(node[1]->coord());
xp[2] = canvas.mesh2xpoint(node[2]->coord());
xp[3] = xp[0];
XDrawLines(gfxinfo.display(), canvas.pixmap(), canvas.gc(), xp, 4,
CoordModeOrigin);
}
void AMTriangle::draw_selected(ImageCanvas &canvas) const {
if(selected()) {
XPoint xp[3];
xp[0] = canvas.mesh2xpoint(node[0]->coord());
xp[1] = canvas.mesh2xpoint(node[1]->coord());
xp[2] = canvas.mesh2xpoint(node[2]->coord());
XFillPolygon(gfxinfo.display(), canvas.pixmap(), canvas.gc(), xp, 3,
Convex, CoordModeOrigin);
}
}
void AMTriangle::fill(ImageCanvas &canvas) const {
XSetForeground(gfxinfo.display(), canvas.gc(), color.pixel);
XPoint xp[3];
xp[0] = canvas.mesh2xpoint(node[0]->coord());
xp[1] = canvas.mesh2xpoint(node[1]->coord());
xp[2] = canvas.mesh2xpoint(node[2]->coord());
XFillPolygon(gfxinfo.display(), canvas.pixmap(), canvas.gc(), xp, 3,
Convex, CoordModeOrigin);
}
void AMTriangle::outline(ImageCanvas &canvas, const Color &outlinecolor) const {
XSetForeground(gfxinfo.display(), canvas.gc(), outlinecolor.pixel);
XSetLineAttributes(gfxinfo.display(), canvas.gc(),
4, LineSolid, CapButt, JoinBevel); // 4 is line width
XPoint xp[4];
xp[0] = canvas.mesh2xpoint(node[0]->coord());
xp[1] = canvas.mesh2xpoint(node[1]->coord());
xp[2] = canvas.mesh2xpoint(node[2]->coord());
xp[3] = xp[0];
XDrawLines(gfxinfo.display(), canvas.pixmap(), canvas.gc(), xp, 4,
CoordModeOrigin);
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
void AMTriangle::writegoof(FILE *file) {
Material *matl = material();
unsigned char flag = 0;
if(matl)
matl->output(file, flag, node[0]->index, node[1]->index, node[2]->index);
else
defaultmaterial->output(file, flag,
node[0]->index, node[1]->index, node[2]->index);
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// elect the most popular material
void AMTriangle::elect_material(Elector &elector) {
if(representative_material_time > mesh->material_recompute_requested &&
representative_material_time > lastmod_time) return;
for(AMTriangleIterator it(*this); !it.end(); ++it) {
const Cell_coordinate pixel = (*this)[it];
elector.vote((*this)[it], intersection(pixel));
}
if(elector.turnout() == 0.0) {
cerr << "There's something wrong!" << endl;
for(AMTriangleIterator it2(*this); !it2.end(); ++it2) {
Cell_coordinate pixel = (*this)[it2];
cerr << pixel << " " << intersection(pixel) << endl;
}
}
elector.tally();
// material_cell = elector.winner();
representative_material = mesh->goof->material[elector.winner()];
++representative_material_time;
}
void AMTriangle::choose_center_material() {
if(representative_material_time > mesh->material_recompute_requested &&
representative_material_time > lastmod_time) return;
// material_cell = closest_pixel();
representative_material = mesh->goof->material[closest_pixel()];
++representative_material_time;
}
Material *AMTriangle::material() {
inherit_material();
return representative_material;
}
Material *AMTriangle::resolve_material() {
Material *mat = material();
if(!mat)
return defaultmaterial;
return mat;
}
// Cell_coordinate AMTriangle::representative_material_cell() {
// inherit_material();
// return material_cell;
// }
void AMTriangle::inherit_material(bool forced) {
if(forced || !inhibit_inheritance || !representative_material) {
if(representative_material_time > mesh->material_recompute_requested &&
representative_material_time > lastmod_time && !forced) return;
switch(int(materialtransfermethod)) {
case MTM_ELECTION:
elect_material(material_elector);
break;
case MTM_CENTERPIXEL:
choose_center_material();
break;
}
}
}
void AMTriangle::set_material(Material *mat, bool inhibit) {
++representative_material_time;
inhibit_inheritance = inhibit;
representative_material = mat;
}
// is this triangle's material different from any of its neighbors?
bool AMTriangle::is_interface() {
Material *my_material = material();
if(!my_material) { // this triangle has no material ...
for(int i=0; i<3; i++)
if(neighbor[i] && neighbor[i]->material())
return 1; // ... but neighbor has a material
return 0; // this tri and all nbrs have no material
}
// this triangle has a material assigned to it
for(int i=0; i<3; i++)
if(neighbor[i])
if(!neighbor[i]->material() || // nbr has no material ...
(neighbor[i]->material() && // ... or nbr has a different material
!(*my_material == *neighbor[i]->material())))
return 1;
return 0;
}
// is this triangle's material different from at least two of its neighbors?
bool AMTriangle::is_double_interface() {
Material *my_material = material();
int nnbrs = 0; // number of neighboring triangles
int ndiff = 0; // number which have different materials
for(int i=0; i<3; i++) {
if(neighbor[i]) {
nnbrs++;
Material *nbr_material = neighbor[i]->material();
if(my_material) {
if(!nbr_material || !(*my_material == *nbr_material))
ndiff++;
}
else { // no material assigned to this triangle
if(nbr_material)
ndiff++;
}
}
}
if(nnbrs == 1 && ndiff == 1) return 1;
if(nnbrs == 2 && ndiff >= 1) return 1;
if(nnbrs == 3 && ndiff >= 2) return 1;
return 0;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// Inherit groups from pixels
void AMTriangle::elect_groups(Elector &elector) {
if(group_cell_time > mesh->groups_recompute_requested) return;
for(AMTriangleIterator it(*this); !it.end(); ++it) {
const Cell_coordinate pixel = (*this)[it];
elector.vote((*this)[it], intersection(pixel));
}
elector.tally();
++group_cell_time;
group_cell = elector.winner();
}
void AMTriangle::copy_groups(const Cell_coordinate &pxl) {
const LinkList<PixelGroup*> &grplist = mesh->goof->pixelgrouplist[pxl];
for(LinkListIterator<PixelGroup*> j = grplist.begin(); !j.end(); ++j) {
// find or create corresponding group in mesh
MeshGroup *g = mesh->find_group(grplist[j]->query_name());
// add this triangle to the group
g->append(this);
}
++meshgroups_time;
}
void AMTriangle::inherit_elected_groups(Elector &elector) {
elect_groups(elector); // computes group_cell
copy_groups(group_cell);
}
void AMTriangle::inherit_center_groups() {
copy_groups(closest_pixel());
}
void AMTriangle::inherit_all_groups() {
// list of all groups contained in any pixel
Vec<MeshGroup*> allgroups(0, BlockSize(100));
for(AMTriangleIterator it(*this); !it.end(); ++it) {
const LinkList<PixelGroup*> &grplist
= mesh->goof->pixelgrouplist[(*this)[it]];
// loop over all groups in pixel
for(LinkListIterator<PixelGroup*> j = grplist.begin(); !j.end(); ++j) {
// see if there are any new groups
int found = 0;
const CharString &grpname = grplist[j]->query_name();
for(int k=0; k<allgroups.capacity() && !found; k++) {
if(grpname == allgroups[k]->query_name())
found = 1;
}
// add new group to allgroups, and also to Mesh::grouplist
if(!found)
allgroups.grow(1, mesh->find_group(grpname));
}
}
// put this triangle in allgroups
for(int n=0; n<allgroups.capacity(); n++)
allgroups[n]->append(this);
++meshgroups_time;
}
void AMTriangle::inherit_groups(bool forced) {
if(forced || !inhibit_groupinheritance) {
if(meshgroups_time > mesh->groups_recompute_requested && !forced)
return;
switch(int(grouptransfermethod)) {
case GTM_ELECTION:
inherit_elected_groups(group_elector);
break;
case GTM_ALLPIXELS:
inherit_all_groups();
break;
case GTM_CENTERPIXEL:
inherit_center_groups();
break;
}
}
}
// static function used by AdaptiveMesh::remove_group() to see if a
// triangle can be removed from a group.
int AMTriangle::inhibittest( AMTriangle *&tri) {
return !tri->inhibit_groupinheritance;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// Add a group to the triangle, but not if it's already there. This
// could be done more efficiently by keeping the list sorted, but is
// it worth the overhead? We don't expect triangles to be in very many
// groups.
void AMTriangle::add_group(MeshGroup *grp) {
const CharString &name = grp->query_name();
for(int i=0; i<meshgroups.capacity(); i++) {
if(meshgroups[i]->query_name() == name)
return;
}
meshgroups.grow(1, grp);
}
void AMTriangle::remove_group(MeshGroup *grp) {
meshgroups.remove(grp);
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// stuff used by MeshGroup
int operator==(const AMTriangle &t1, const AMTriangle &t2) {
if(*t1.node[0] != *t2.node[0]) return 0;
if(*t1.node[1] != *t2.node[1]) return 0;
if(*t1.node[2] != *t2.node[2]) return 0;
return 1;
}
int operator<(const AMTriangle &t1, const AMTriangle &t2) {
if(*t1.node[0] < *t2.node[0]) return 1;
if(*t1.node[0] > *t2.node[0]) return 0;
if(*t1.node[1] < *t2.node[1]) return 1;
if(*t1.node[1] > *t2.node[1]) return 0;
if(*t1.node[2] < *t2.node[2]) return 1;
return 0;
}
template <>
void MeshGroup::append(AMTriangle* const &tri) {
append_base(tri);
tri->add_group(this);
}
template <>
void MeshGroup::append(const Vec<AMTriangle*> &list) {
append_base(list);
for(int i=0; i<list.capacity(); i++)
list[i]->add_group(this);
}
template <>
void MeshGroup::clear() {
for(int i=0; i<members.capacity(); i++)
members[i]->meshgroups.remove(this);
clear_base();
}
template <>
void MeshGroup::remove(AMTriangle * const &tri) {
tri->remove_group(this);
remove_base(tri);
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
Vec<Cell_coordinate> AMTriangle::get_pixels() const {
Vec<Cell_coordinate> list(0, BlockSize(10));
for(AMTriangleIterator it(*this); !it.end(); ++it)
list.grow(1, (*this)[it]);
return list;
}
int AMTriangle::npixels() const {
int n = 0;
for(AMTriangleIterator it(*this); !it.end(); ++it)
n++;
return n;
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
AMTriangle *AMTriangle::copy(AdaptiveMesh *newmesh) const {
AMTriangle *newtri = new AMTriangle;
int i;
for(i=0; i<2; i++) {
if(child[i]) {
newtri->child[i] = child[i]->copy(newmesh);
newtri->child[i]->parent = newtri;
}
else {
newtri->child[0] = 0;
newtri->child[1] = 0;
}
}
newtri->mesh = newmesh;
newtri->generation = generation;
newtri->index = index;
newtri->swapped = swapped;
newtri->status_ = 0;
if(selected())
newtri->select();
newtri->parent = 0; // will be set by parent
newtri->representative_material = representative_material;
// newtri->material_cell = material_cell;
newtri->inhibit_inheritance = inhibit_inheritance;
newtri->representative_material_time = representative_material_time;
newtri->group_cell = group_cell;
newtri->group_cell_time = group_cell_time;
newtri->meshgroups_time = meshgroups_time;
newtri->inhibit_groupinheritance = inhibit_groupinheritance;
newtri->lastE = lastE;
newtri->currentE = currentE;
newtri->Ecalculated = Ecalculated;
newtri->lastEcalculated = lastEcalculated;
for(i=0; i<3; i++)
newtri->node[i] = newmesh->nodes[node[i]->index];
if(!child[0])
// Add triangle to list in nodes iff it has no children.
for(i=0; i<3; i++)
newtri->node[i]->add_triangle(newtri);
// Can't set neighbor pointers here, because neighbor may not have
// been copied yet.
for(i=0; i<3; i++)
newtri->neighbor[i] = 0;
for(i=0; i<meshgroups.capacity(); i++) {
MeshGroup *grp = newmesh->get_group(meshgroups[i]->query_name());
// newtri->meshgroups.grow(1, grp);
if(grp) // don't worry about unnamed groups
grp->append(newtri); // adds grp to newtri->meshgroups too
}
return newtri;
}
void AMTriangle::add_neighbor(int n0, int n1, AMTriangle *nbr) {
// n0 and n1 are node indices. nbr is the neighbor opposite the OTHER node.
for(int i=0; i<3; i++) {
if(node[i]->index == n0) {
if(node[(i+1)%3]->index == n1) {
#ifdef DEBUG
if(neighbor[(i+2)%3] != 0)
cerr << "Error 1 in AMTriangle::add_neighbor()! "
<< neighbor[(i+2)%3]
<< endl;
#endif
neighbor[(i+2)%3] = nbr;
}
#ifndef DEBUG
else {
#else
else if(node[(i+2)%3]->index == n1) {
if(neighbor[(i+1)%3] != 0)
cerr << "Error 2 in AMTriangle::add_neighbor()! "
<< neighbor[(i+1)%3]
<< endl;
#endif
neighbor[(i+1)%3] = nbr;
}
#ifdef DEBUG
else
cerr << "Error 3 in AMTriangle::add_neighbor()! " << endl;
#endif
return;
}
}
cerr << "Error in AMTriangle::add_neighbor()!" << endl;
exit(1);
}
//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//=\\=//
// The Sun compiler seems to want these declared explicitly...
#if defined(sun) && !defined(__GNUG__)
template class LinkList<Cell_coordinate>::LinkListNode;
template class LinkList<AMTriangle*>::LinkListNode;
#endif
syntax highlighted by Code2HTML, v. 0.9.1