gridsub2.c

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

void rectblk::initsides(void){

    int i;

    // change rectblk to this style

    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;
    }
}


// destructor: cleanup

rectblk::~rectblk(){

    delete [] exteriornode;
    delete [] exteriorside;
}


// assignment: cleanup and copy

rectblk& rectblk::operator = (rectblk& 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;
      length1 = a.length1;
      length2 = a.length2;
      div1    = a.div1;
      div2    = a.div2;
      angle   = a.angle;
      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 rectblk::draw(HDC hdc){

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

void rectblk::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 rectblk :: highlight(HDC hdc){

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

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

    outline(hdc, nDrawingMode, hpen, length1, length2, dpoint, 0.0, 1.0);
}

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

    outline(hdc, nDrawingMode, hpen, length1, length2, center, sintheta, costheta);
}


// outline resized object
void rectblk::outline(HDC hdc, short nDrawingMode, HPEN hpen, dPOINT corner, int hookedcorner){

    dPOINT newloc;
    double newl1, newl2;

    eval(corner, &newloc, &newl1, &newl2, hookedcorner);
    outline(hdc, nDrawingMode, hpen, newl1, newl2, newloc, 0.0, 1.0);
}


void rectblk::outline(HDC hdc, short nDrawingMode, HPEN hpen, double factor, dPOINT loc){

    double newl1, newl2;

    newl1 = factor * getl1();
    newl2 = factor * getl2();
    outline(hdc, nDrawingMode, hpen, newl1, newl2, loc, 0.0, 1.0);
}


// outline kernel
void rectblk::outline(HDC hdc, short nDrawingMode, HPEN hpen, double width, double height, dPOINT center, double sintheta, double costheta){

//  center is the center about which the object is rotated.
//  theta is the ADDITIONAL angle to rotate.
//  short  nDrawingMode = R2_NOTXORPEN;
    dPOINT p1, p2, p3, p4;
    double x1, y1, x2, y2, x3, y3, x4, y4;
    double sinalpha, cosalpha, sinbeta, cosbeta;

//  alpha = beta + theta

    sinbeta = sin(angle * PI / 180.0);
    cosbeta = cos(angle * PI / 180.0);

    sinalpha = sinbeta * costheta + cosbeta * sintheta;
    cosalpha = cosbeta * costheta - sinbeta * sintheta;

    p1.x = 0.0;
    p1.y = 0.0;
    x1 = p1.x + center.x;
    y1 = p1.y + center.y;

    p2.x = width;
    p2.y = 0.0;
    rotatepoint(&p2, sinalpha, cosalpha);
    x2 = p2.x + center.x;
    y2 = p2.y + center.y;

    p3.x = + width;
    p3.y = + height;
    rotatepoint(&p3, sinalpha, cosalpha);
    x3 = p3.x + center.x;
    y3 = p3.y + center.y;

    p4.x = 0.0;
    p4.y = + height;
    rotatepoint(&p4, sinalpha, cosalpha);
    x4 = p4.x + center.x;
    y4 = p4.y + center.y;

    // 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 * rectblk::copy(void){

    int i;
    rectblk *newblock;

    newblock = new rectblk();

    (*newblock) = (*this);

    return newblock;
}


void rectblk::resize(double l1, double l2){

    length1 = l1;
    length2 = l2;
    snaptogrid();
    generatenodes();

}

void rectblk::resize(double factor){

    double l1, l2;

    l1 = factor * getl1();
    l2 = factor * getl2();
    resize(l1, l2);
}

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

    dPOINT newloc;
    double newl1, newl2;

    eval(corner, &newloc, &newl1, &newl2, hookedcorner);
    block::move(newloc.x, newloc.y);
    resize(newl1, newl2);
}


void rectblk::refine(){

}


void rectblk::reconfig(double ep, double sgm, double m, double m_sgm,
		 double x, double y, double l1, double l2,
		 int d1, int d2, double ang){

    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;
    length1  = l1;
    length2  = l2;
    angle    = ang;
}



void rectblk::dialog(void){

    rectblockselected = this;
    SendMessage(hwndMain, WM_COMMAND, IDM_RECTANGLE, 0L);

}


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

    dPOINT p;
    double sintheta, costheta;

    sintheta = sin(angle * PI / 180.0);
    costheta = cos(angle * PI / 180.0);


    switch(corner){
      case 1:
	p.x = length1;
	p.y = length2;
      break;

      case 2:
	p.x = 0.0;
	p.y = length2;
      break;

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

      case 4:
	p.x = length1;
	p.y = 0.0;
      break;

      default:
	dpoint->x = 0.0;
	dpoint->y = 0.0;
    }
    rotatepoint(&p, sintheta, costheta);
    dpoint->x = xloc + p.x;
    dpoint->y = yloc + p.y;

}


void rectblk::eval(dPOINT dpoint, dPOINT *newloc, double *newwidth, double *newheight, int corner){

    double sintheta, costheta;
    dPOINT p;

    sintheta = sin(-angle * PI / 180.0);
    costheta = cos(-angle * PI / 180.0);

    p.x = dpoint.x - xloc;
    p.y = dpoint.y - yloc;
    rotatepoint(&p, sintheta, costheta);

    switch(corner){
      case 1:
	*newwidth  = p.x;
	if(*newwidth < 0.0) *newwidth = 0;
	*newheight = p.y;
	if(*newheight < 0.0) *newheight = 0;
	newloc->x = 0.0;
	newloc->y = 0.0;
      break;

      case 2:
	*newwidth  = length1 - p.x;
	if(*newwidth < 0.0) *newwidth = 0;
	*newheight = p.y;
	if(*newheight < 0.0) *newheight = 0;
	newloc->x = p.x;
	newloc->y = 0.0;
      break;

      case 3:
	*newwidth  = length1 - p.x;
	if(*newwidth < 0.0) *newwidth = 0;
	*newheight = length2 - p.y;
	if(*newheight < 0.0) *newheight = 0;
	newloc->x = p.x;
	newloc->y = p.y;
      break;

      case 4:
	*newwidth  = p.x;
	if(*newwidth < 0.0) *newwidth = 0;
	*newheight = length2 - p.y;
	if(*newheight < 0.0) *newheight = 0;
	newloc->x = 0.0;
	newloc->y = p.y;
      break;
    }
    rotatepoint(newloc, -sintheta, costheta);
    newloc->x += xloc;
    newloc->y += yloc;
}

double rectblk::allowablesnap(void){

    double snap;

    snap = 0.5 * length1 / div1;
    if(0.5 * length2 / div2 < snap) snap = 0.5 * length2 / div2;
    return snap;
}


void rectblk::snaptogrid(void){

    if(option.checkposition == YES){
      xloc = (double)(option.position * floor(0.5 + xloc/option.position));
      yloc = (double)(option.position * floor(0.5 + yloc/option.position));
      length1 = (double)(option.position * floor(0.5 + length1/option.position));
      length2 = (double)(option.position * floor(0.5 + length2/option.position));
      if(length1 < option.position/2.0){
	length1 = option.position;
      }
      if(length2 <= option.position/2.0){
	length2 = option.position;
      }
    }
    if(option.checkrotation == YES){
      angle = (double)(option.rotation * floor(0.5+angle/option.rotation));
    }
}



COUNT rectblk::gencellinfo(FILE *stream, COUNT count){

    unsigned int i, j, k; // elm;
    unsigned int side1, side2, side3, side4;
    COUNT n;

    n = count;
    for(j = 1; j <= div2; j++)
      for (i = 1; i <= div1; i++){
	fprintf(stream, "%u %u %lg %lg %lg %lg ",
		n.cell, RECTCODE, epsilon, sigma, mu, sigma_mu);

	if(j == 1){    //side 1
	  side1 = exteriorside[i-1].index;
	}
	else{
	  side1 = i + (j-2) * div1 + n.side - 1;
	}

	if(i == div1){  // side 2
	  side2 = exteriorside[div1 + j - 1].index;
	}
	else{
	  side2 = j + (i-1) * div2 + (div2 - 1) * div1 + n.side - 1;
	}

	if(j == div2){  // side 3
	  side3 = exteriorside[2 * div1 + div2 - i].index;
	}
	else{
	  side3 = i + (j - 1) * div1 + n.side - 1;
	}

	if(i == 1){    // side 4
	  side4 = exteriorside[2 * (div1 + div2) - j].index;
	}
	else{
	  side4 = j + (i-2) * div2 + (div2-1) * div1 + n.side - 1;
	}
	if(EOF == fprintf(stream, "%u %u %u %u\n", side1, side2, side3, side4)){
	  n.node = 0;
	  n.side = 0;
	  return n;
	}

	(n.cell)++;
      }

    n.side += (div1 - 1) * div2 + (div2 - 1) * div1;
    n.node += (div1 - 1) * (div2 - 1);
    return n;
}


COUNT rectblk::geninteriorsideinfo(FILE *stream, COUNT count){

    unsigned int i, j, k; // elm;
    unsigned int sindex, sn, en;
    COUNT n;

    n = count;
    for(j = 1; j < div2; j++)
      for (i = 1; i <= div1; i++){
	sindex = i + (j - 1) * div1 + n.side - 1;
	if(i == 1)
	  sn = exteriornode[2 * (div1 + div2) - j].index;
	else
	  sn = (i-1) + (j-1) * (div1-1) + n.node - 1;
	if(i == div1)
	  en = exteriornode[div1 + j].index;
	else
	  en = i + (j-1) * (div1-1) + n.node - 1;

	if(EOF == fprintf(stream, "%u %u %u %u\n", sindex, 0, sn, en)){
	  n.node = 0;
	  n.side = 0;
	  return n;
	}
      }

    for (i = 1; i < div1; i++)
      for(j = 1; j <= div2; j++){
	sindex = j + (i-1) * div2 + (div2 - 1) * div1 + n.side - 1;
	if(j == 1)
	  sn = exteriornode[i].index;
	else
	  sn = i + (j-2) * (div1-1) + n.node - 1;

	if(j == div2)
	  en = exteriornode[2 * div1 + div2 - i].index;