gridsub2.c

有限元学习研究用源代码(老外的),供科研人员参考

      p.x = (1.0 - (double)i / div2) * p2.x + ((double)i / div2) * p3.x;
      p.y = (1.0 - (double)i / div2) * p2.y + ((double)i / div2) * p3.y;
      exteriornode[2 * div1 + div2 + i].x = p.x;
      exteriornode[2 * div1 + div2 + i].y = p.y;
    }
}

void quadblk::initnodes(void){

    int i;
    for(i = 0; i < 2 * (div1 + div2); i++){
      exteriornode[i].cellid = id;
      exteriornode[i].node = i;
    }
}


void quadblk::initsides(void){

    int i;

    for(i = 0; i < div1; i++){
      exteriorside[i].cellid = id;
      exteriorside[i].side = i;
    }

    for(i = 0; i < div2; i++){
      exteriorside[div1 + i].cellid = id;
      exteriorside[div1 + i].side = div1 + i;
    }

    for(i = 0; i < div1; i++){
      exteriorside[div1 + div2 + i].cellid = id;
      exteriorside[div1 + div2 + i].side = div1 + div2 + i;
    }

    for(i = 0; i < div2; i++){
      exteriorside[2 * div1 + div2 + i].cellid = id;
      exteriorside[2 * div1 + div2 + i].side = 2 * div1 + div2 + i;
    }
}


void quadblk::move(dPOINT center, double sintheta, double costheta){

//  center is the center at which the object is to be outlined.
//  theta is the ADDITIONAL angle to rotate.
    dPOINT p;

    p.x = cornerx1;
    p.y = cornery1;
    rotatepoint(&p, sintheta, costheta);
    cornerx1 = p.x;
    cornery1 = p.y;

    p.x = cornerx2;
    p.y = cornery2;
    rotatepoint(&p, sintheta, costheta);
    cornerx2 = p.x;
    cornery2 = p.y;

    p.x = cornerx3;
    p.y = cornery3;
    rotatepoint(&p, sintheta, costheta);
    cornerx3 = p.x;
    cornery3 = p.y;

    xloc = center.x;
    yloc = center.y;
    snaptogrid();
    generatenodes();
}


// destructor: cleanup

quadblk::~quadblk(){

    delete [] exteriornode;
    delete [] exteriorside;
}


// assignment: cleanup and copy

quadblk& quadblk::operator = (const quadblk& a)
{
    int i, elmt;

    if(this != &a){
      delete [] exteriornode;
      delete [] exteriorside;
      epsilon  = a.epsilon;
      sigma    = a.sigma;
      mu       = a.mu;
      sigma_mu = a.sigma_mu;
      xloc    = a.xloc;
      yloc    = a.yloc;
      cornerx1 = a.cornerx1;
      cornery1 = a.cornery1;
      cornerx2 = a.cornerx2;
      cornery2 = a.cornery2;
      cornerx3 = a.cornerx3;
      cornery3 = a.cornery3;
      div1    = a.div1;
      div2    = a.div2;
      id =  a.id;
      elmt = exteriorelmtcount();
      exteriornode = new surnode[elmt];
      exteriorside = new surside[elmt];

      for(i = 0; i < elmt; i++){
	exteriornode[i] = a.exteriornode[i];
	exteriorside[i] = a.exteriorside[i];
      }

    }
    return *this;
}

void quadblk::draw(HDC hdc){

    if(outsidewindow()) return;
    drawinteriorcell(hdc, hBrushNormal);
    drawsidesandnodes(hdc);
}

void quadblk::drawinteriorcell(HDC hdc, HBRUSH hBrush){

    int i;
    double xs, ys, xe, ye, xf, yf;
    dPOINT dpt;
    POINT pt[5];

    vertex(&dpt, 1);
    drawingtolp(&pt[0], dpt);
    drawingtolp(&pt[4], dpt);
    vertex(&dpt, 2);
    drawingtolp(&pt[1], dpt);
    vertex(&dpt, 3);
    drawingtolp(&pt[2], dpt);
    vertex(&dpt, 4);
    drawingtolp(&pt[3], dpt);

    SetPolyFillMode(hdc, WINDING);
    SelectObject(hdc, hPenSolidYellow);
    SelectObject(hdc, hBrush);
    SetROP2 (hdc, R2_NOTXORPEN);
    Polygon(hdc, pt, 5);

    SelectObject(hdc, hPenSolidBlack);
    SetROP2 (hdc, R2_COPYPEN) ;

    if((div1 != 1) || (div2 != 1)){
      for(i = 1; i < div1; i++){
	drawline(hdc, exteriornode[i].x, exteriornode[i].y,
	  exteriornode[2*div1 + div2 - i].x,
	  exteriornode[2*div1 + div2 - i].y);
      }

      for(i = div1+1; i < div1 + div2; i++){
	drawline(hdc, exteriornode[i].x, exteriornode[i].y,
	  exteriornode[3*div1 + 2*div2 - i].x,
	  exteriornode[3*div1 + 2*div2 - i].y);
      }
    }
}


void quadblk :: highlight(HDC hdc){

    if(outsidewindow()) return;
    drawinteriorcell(hdc, hBrushHighlight);
    drawsidesandnodes(hdc);
}

// 1. outline moved object
void quadblk::outline(HDC hdc, short nDrawingMode, HPEN hpen, dPOINT  dpoint){

    dPOINT newcorner = {0.0, 0.0};
    outline(hdc, nDrawingMode, hpen, 1.0, dpoint, newcorner, 0, 0.0, 1.0);
}

// 2. outline rotated object
void quadblk::outline(HDC hdc, short nDrawingMode, HPEN hpen, dPOINT center, double sintheta, double costheta){

    dPOINT newcorner = {0.0, 0.0};
    outline(hdc, nDrawingMode, hpen, 1.0, center, newcorner, 0, sintheta, costheta);
}


// 3. outline reshaped object
void quadblk::outline(HDC hdc, short nDrawingMode, HPEN hpen, dPOINT corner, int hookedcorner){

    dPOINT center;
    getreference(&center);
    outline(hdc, nDrawingMode, hpen, 1.0, center, corner, hookedcorner, 0.0, 1.0);
}

// 4. outline rescaled object
void quadblk::outline(HDC hdc, short nDrawingMode, HPEN hpen, double scale, dPOINT center){

    dPOINT dummycorner = {0.0, 0.0};
    outline(hdc, nDrawingMode, hpen, scale, center, dummycorner, 0, 0.0, 1.0);
}


// 5. outline kernel
void quadblk::outline(HDC hdc, short nDrawingMode, HPEN hpen, double scale, dPOINT center, dPOINT newcorner, int hookedcorner, double sintheta, double costheta){

//  scale  is the resize factor
//  newcorner is the new position the hooked corener to be moved to.
//  center is the center at which the object is to be outlined.
//  theta is the ADDITIONAL angle to rotate.
//  short  nDrawingMode = R2_NOTXORPEN;
    dPOINT p, p1, p2, p3, p4;
    double x1, y1, x2, y2, x3, y3, x4, y4;

    vertex(&p, 1);
    p1.x = p.x - xloc;
    p1.y = p.y - yloc;
    rotatepoint(&p1, sintheta, costheta);
    x1 = scale * p1.x + center.x;
    y1 = scale * p1.y + center.y;


    vertex(&p, 2);
    p2.x = p.x - xloc;
    p2.y = p.y - yloc;
    rotatepoint(&p2, sintheta, costheta);
    x2 = scale * p2.x + center.x;
    y2 = scale * p2.y + center.y;

    vertex(&p, 3);
    p3.x = p.x - xloc;
    p3.y = p.y - yloc;
    rotatepoint(&p3, sintheta, costheta);
    x3 = scale * p3.x + center.x;
    y3 = scale * p3.y + center.y;

    vertex(&p, 4);
    p4.x = p.x - xloc;
    p4.y = p.y - yloc;
    rotatepoint(&p4, sintheta, costheta);
    x4 = scale * p4.x + center.x;
    y4 = scale * p4.y + center.y;

    switch(hookedcorner){

       case 1:
	 x1 = newcorner.x;
	 y1 = newcorner.y;
       break;

       case 2:
	 x2 = newcorner.x;
	 y2 = newcorner.y;
       break;

       case 3:
	 x3 = newcorner.x;
	 y3 = newcorner.y;
       break;

       case 4:
	 x4 = newcorner.x;
	 y4 = newcorner.y;
       break;
    }

    // Draw the outline with a black pen.

    SelectObject (hdc, hpen) ;
    SetROP2 (hdc, nDrawingMode) ;

    drawline(hdc, x1, y1, x2, y2);
    drawline(hdc, x2, y2, x3, y3);
    drawline(hdc, x3, y3, x4, y4);
    drawline(hdc, x4, y4, x1, y1);
}

block * quadblk::copy(void){

    int i;
    quadblk *newblock;

    newblock = new quadblk();

    *newblock = *this;

    return newblock;
}


void quadblk::resize(double factor){

    double l1, l2, l3;

    cornerx1 *= factor;
    cornery1 *= factor;
    cornerx2 *= factor;
    cornery2 *= factor;
    cornerx3 *= factor;
    cornery3 *= factor;
    snaptogrid();
    generatenodes();
}

void quadblk::resize(dPOINT corner, int hookedcorner){

    dPOINT corner1, corner2, corner3, corner4;
    dPOINT v1, v2;

    switch(hookedcorner){
      case 1:
	cornerx2 = corner.x - xloc;
	cornery2 = corner.y - yloc;
      break;

      case 2:
	cornerx3 = corner.x - xloc;
	cornery3 = corner.y - yloc;
      break;

      case 3:
	vertex(&corner1, 1);
	vertex(&corner2, 2);
	vertex(&corner4, 4);
	xloc = corner.x;
	yloc = corner.y;
	cornerx1 = corner4.x - xloc;
	cornery1 = corner4.y - yloc;
	cornerx2 = corner1.x - xloc;
	cornery2 = corner1.y - yloc;
	cornerx3 = corner2.x - xloc;
	cornery3 = corner2.y - yloc;
      break;

      case 4:
	cornerx1 = corner.x - xloc;
	cornery1 = corner.y - yloc;
      break;
    }
    ordervertex();
    snaptogrid();
    generatenodes();
}


void quadblk::refine(){

}


void quadblk::reconfig(double ep, double sgm, double m, double m_sgm,
		 double x, double y, double x1, double y1,
		 double x2, double y2, double x3, double y3,
		 int d1, int d2){

    if((div1 != d1) || (div2 != d2)){
      delete [] exteriornode;
      delete [] exteriorside;
      div1     = d1;
      div2     = d2;
      exteriornode = new surnode[exteriorelmtcount()];
      exteriorside = new surside[exteriorelmtcount()];
    }
    epsilon  = ep;
    sigma    = sgm;
    mu       = m;
    sigma_mu = m_sgm;
    xloc     = x;
    yloc     = y;
    cornerx1 = x1;
    cornery1 = y1;
    cornerx2 = x2;
    cornery2 = y2;
    cornerx3 = x3;
    cornery3 = y3;
}



void quadblk::dialog(void){

    quadblockselected = this;
    SendMessage(hwndMain, WM_COMMAND, IDM_QUADRILATERAL, 0L);

}


void quadblk::vertex(dPOINT *dpoint, int corner){

    dPOINT p;


    switch(corner){
      case 1:
	p.x = cornerx2;
	p.y = cornery2;
      break;

      case 2:
	p.x = cornerx3;
	p.y = cornery3;
      break;

      case 3:
	p.x = 0.0;
	p.y = 0.0;
      break;

      case 4:
	p.x = cornerx1;
	p.y = cornery1;
      break;

      default:
	dpoint->x = 0.0;
	dpoint->y = 0.0;
    }
    dpoint->x = xloc + p.x;
    dpoint->y = yloc + p.y;
}

double quadblk::allowablesnap(void){

    double snap, snap1, snap2, snap3, snap4;
    double snap5, snap6, snap7, snap8;
    dPOINT c1, c2, c3, c4, cs, ce;

    vertex(&c1, 1);
    vertex(&c2, 2);
    vertex(&c3, 3);
    vertex(&c4, 4);
    snap1 = 0.5 * length(c1, c2) / div1;
    snap2 = 0.5 * length(c2, c3) / div2;
    snap3 = 0.5 * length(c3, c4) / div1;
    snap4 = 0.5 * length(c4, c1) / div2;
    cs.x = (c2.x - c3.x) / div2;
    cs.y = (c2.y - c3.y) / div2;
    ce.x = (c4.x - c3.x) / div1;
    ce.y = (c4.y - c3.y) / div1;
    snap5 = 0.5 * length(cs, ce);
    cs.x = (c1.x - c4.x) / div2;
    cs.y = (c1.y - c4.y) / div2;
    ce.x = (c3.x - c4.x) / div1;
    ce.y = (c3.y - c4.y) / div1;
    snap6 = 0.5 * length(cs, ce);
    cs.x = (c4.x - c1.x) / div2;
    cs.y = (c4.y - c1.y) / div2;
    ce.x = (c2.x - c1.x) / div1;
    ce.y = (c2.y - c1.y) / div1;
    snap7 = 0.5 * length(cs, ce);
    cs.x = (c3.x - c2.x) / div2;
    cs.y = (c3.y - c2.y) / div2;
    ce.x = (c1.x - c2.x) / div1;
    ce.y = (c1.y - c2.y) / div1;
    snap8 = 0.5 * length(cs, ce);


    snap = snap1 < snap2 ? snap1 : snap2;
    snap = snap < snap3 ? snap : snap3;
    snap = snap < snap4 ? snap : snap4;
    snap = snap < snap5 ? snap : snap5;
    snap = snap < snap6 ? snap : snap6;
    snap = snap < snap7 ? snap : snap7;
    snap = snap < snap8 ? snap : snap8;

    // remember to take into account incline angle
    return snap;
}


void quadblk::snaptogrid(void){

    double x1, y1, x2, y2, x3, y3;

    if(option.checkposition == NO) return;


    xloc = (double)(option.position * floor(0.5 + xloc/option.position));
    yloc = (double)(option.position * floor(0.5 + yloc/option.position));

    x1 = cornerx1 + xloc;
    y1 = cornery1 + yloc;

    x2 = cornerx2 + xloc;
    y2 = cornery2 + yloc;

    x3 = cornerx3 + xloc;
    y3 = cornery3 + yloc;


    x1 = (double)(option.position * floor(0.5 + x1/option.position));
    y1 = (double)(option.position * floor(0.5 + y1/option.position));

    x2 = (double)(option.position * floor(0.5 + x2/option.position));
    y2 = (double)(option.position * floor(0.5 + y2/option.position));

    x3 = (double)(option.position * floor(0.5 + x3/option.position));
    y3 = (double)(option.position * floor(0.5 + y3/option.position));


    cornerx1 =  x1 - xloc;
    cornery1 =  y1 - yloc;
    cornerx2 =  x2 - xloc;
    cornery2 =  y2 - yloc;
    cornerx3 =  x3 - xloc;
    cornery3 =  y3 - yloc;

    // Make sure the block is not zero-size.
    // it must has at least a "unit" size.
    if((fabs(cornerx1) < option.