gridsub2.c

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

	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;
	}
      }
    n.side += (div1 - 1) * div2 + (div2 - 1) * div1;
    n.node += (div1 - 1) * (div2 - 1);
    return n;
}


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

    unsigned int    i, j, n;
    unsigned long deno;
    dPOINT p1, p2, p3, p4;
    COUNT  c;
    double x, y;

    deno = long(div1) * long(div2);
    vertex(&p1, 1);
    vertex(&p2, 2);
    vertex(&p3, 3);
    vertex(&p4, 4);
    n = count.node;
    for(j = 1; j < div2; j++)
      for(i = 1; i < div1; i++){
	x = ((p3.x * (div1 - i) + p4.x * i) * (div2 - j) +
	     (p2.x * (div1 - i) + p1.x * i) * j) / deno; 
	y = ((p3.y * (div1 - i) + p4.y * i) * (div2 - j) +
	     (p2.y * (div1 - i) + p1.y * i) * j) / deno; 
	fprintf(stream, "%u %u %lg %lg\n", n, 0, x, y);
	n++;
      }
    c.node = count.node + long(div1 - 1) * long(div2 - 1);
    c.side = count.side;
    return c;
}

int rectblk::savegdfinfo(FILE *stream){

   unsigned int i, j,  k, n;

   n = exteriorelmtcount();

   if(fprintf(stream, "%u %u\n", RECTCODE, id) < 0) return FALSE;
   if(fprintf(stream, "%lg %lg %lg %lg %lg %lg\n",
     epsilon, sigma, mu, sigma_mu, xloc, yloc) < 0) return FALSE;
   if(fprintf(stream, "%lg %lg %lg %u %u\n",
     length1, length2, angle, div1, div2) < 0) return FALSE;
   for(i = 0; i < n; i++){
     if(fprintf(stream, "%u %u %u %u %u %u %u\n",
       i, exteriornode[i].node, exteriornode[i].cellid, exteriornode[i].property & PROPERTYMASK,
	  exteriorside[i].side, exteriorside[i].cellid, exteriorside[i].property & PROPERTYMASK)
	< 0) return FALSE;
   }
   if(fprintf(stream, "\n") < 0) return FALSE;
   return TRUE;
}

int rectblk::readgdfinfo(FILE *stream){

   unsigned int i, n;

   fscanf(stream, "%lf %lf %lf %lf %lf %lf\n",
     &epsilon, &sigma, &mu, &sigma_mu, &xloc, &yloc);
   fscanf(stream, "%lf %lf %lf %u %u\n",
     &length1, &length2, &angle, &div1, &div2);

   n = exteriorelmtcount();

   exteriornode = new surnode[n];
   exteriorside = new surside[n];

   generatenodes();

   for(i = 0; i < n; i++){
     fscanf(stream, "%*u %u %u %u %u %u %u\n",
       &exteriornode[i].node, &exteriornode[i].cellid, &exteriornode[i].property,
       &exteriorside[i].side, &exteriorside[i].cellid, &exteriorside[i].property);
   }
   return TRUE;
}


dPOINT rectblk::snaptonode(dPOINT p){

    dPOINT p1, p2, p3, p4, rp;
    double c1, c2, c3, c4, c5, c6, c7, c8;
    double b, c;
    int    i, j;
    unsigned long deno;

    vertex(&p1, 1);
    vertex(&p2, 2);
    vertex(&p3, 3);
    vertex(&p4, 4);

    deno = long(div1) * long(div2);

    c1 = p1.y - p2.y + p3.y - p4.y;
    c2 = (p4.y - p3.y) * div2;
    c3 = (p2.x - p3.x) * div1;
    c4 = (p3.x - p.x) * deno;

    c5 = p1.x - p2.x + p3.x - p4.x;
    c6 = (p4.x - p3.x) * div2;
    c7 = (p2.y - p3.y) * div1;
    c8 = (p3.y - p.y) * deno;

    b = c1 * c4 + c2 * c3 - c5 * c8 - c6 * c7;
    c = c2 * c4 - c6 * c8;

    j = int(0.5 - c / b);

    c1 = p1.y - p2.y + p3.y - p4.y;
    c2 = (p2.y - p3.y) * div1;
    c3 = (p4.x - p3.x) * div2;
    c4 = (p3.x - p.x) * deno;

    c5 = p1.x - p2.x + p3.x - p4.x;
    c6 = (p2.x - p3.x) * div1;
    c7 = (p4.y - p3.y) * div2;
    c8 = (p3.y - p.y) * deno;

    b = c1 * c4 + c2 * c3 - c5 * c8 - c6 * c7;
    c = c2 * c4 - c6 * c8;

    i = int(0.5 - c / b);

    rp.x = ((p3.x * (div1 - i) + p4.x * i) * (div2 - j) +
	    (p2.x * (div1 - i) + p1.x * i) * j) / deno; 
    rp.y = ((p3.y * (div1 - i) + p4.y * i) * (div2 - j) +
	    (p2.y * (div1 - i) + p1.y * i) * j) / deno; 

    return rp;
}

int rectblk::snaptoCSN(dPOINT p, unsigned int *ii, unsigned int *jj){

    dPOINT p1, p2, p3, p4, rp;
    double c1, c2, c3, c4, c5, c6, c7, c8;
    double b, c;
    int snaptype;
    unsigned int i, j; 
    unsigned long deno;

    vertex(&p1, 1);
    vertex(&p2, 2);
    vertex(&p3, 3);
    vertex(&p4, 4);

    deno = 4 * long(div1) * long(div2);

    c1 = p1.y - p2.y + p3.y - p4.y;
    c2 = (p4.y - p3.y) * 2.0 * div2;
    c3 = (p2.x - p3.x) * 2.0 * div1;
    c4 = (p3.x - p.x) * deno;

    c5 = p1.x - p2.x + p3.x - p4.x;
    c6 = (p4.x - p3.x) * 2.0 * div2;
    c7 = (p2.y - p3.y) * 2.0 * div1;
    c8 = (p3.y -  p.y) * deno;

    b = c1 * c4 + c2 * c3 - c5 * c8 - c6 * c7;
    c = c2 * c4 - c6 * c8;

    j = int(0.5 - c / b);

    c1 = p1.y - p2.y + p3.y - p4.y;
    c2 = (p2.y - p3.y) * 2.0 * div1;
    c3 = (p4.x - p3.x) * 2.0 * div2;
    c4 = (p3.x - p.x) * deno;

    c5 = p1.x - p2.x + p3.x - p4.x;
    c6 = (p2.x - p3.x) * 2.0 * div1;
    c7 = (p4.y - p3.y) * 2.0 * div2;
    c8 = (p3.y - p.y) * deno;

    b = c1 * c4 + c2 * c3 - c5 * c8 - c6 * c7;
    c = c2 * c4 - c6 * c8;


    i = int(0.5 - c / b);

    if((i % 2 == 0) && (j %2 == 0)){
      snaptype = NODESNAPED;
      *ii = i / 2;
      *jj = j / 2;
    }
    else if ((i % 2 == 1) && (j % 2 == 1)){
      snaptype = CELLSNAPED;
      *ii = i / 2 + 1;
      *jj = j / 2 + 1;
    }
    else{
      snaptype = SIDESNAPED;
      *ii = i; 
      *jj = j;
    }
    return snaptype;
}

CELLINFO rectblk::cellenquiry(unsigned int i, unsigned int j){

    CELLINFO cellinfo;

    cellinfo.index    = (j-1) * div1 + i;
    cellinfo.epsilon  = epsilon;
    cellinfo.mu       = mu;
    cellinfo.sigma    = sigma;
    cellinfo.sigma_mu = sigma_mu;

    return cellinfo;
}

SIDEINFO rectblk::sideenquiry(unsigned int i, unsigned int j){

    SIDEINFO sideinfo;

    if (i % 2 == 0){
      i = i / 2;
      j = j / 2 + 1;
      if(i == 0){
	sideinfo.where = EXTERIOR;
	sideinfo.index = exteriorside[2*(div1+div2) - j].index;
	sideinfo.property = exteriorside[2*(div1+div2) - j].property;
      }
      else if (i == div1){
	sideinfo.where = EXTERIOR;
	sideinfo.index = exteriorside[div1-1+j].index;
	sideinfo.property = exteriorside[div1-1+j].property;
      }
      else{
	sideinfo.where = INTERIOR;
	sideinfo.index = div1 * (div2-1) +
			 (i-1) * (div2) + j;
	sideinfo.property = INTERIOR;
      }
    }
    else{ // j % 2 == 0
      i = i / 2 + 1;
      j = j / 2;
      if(j == 0){
	sideinfo.where = EXTERIOR;
	sideinfo.index = exteriorside[i-1].index;
	sideinfo.property = exteriorside[i-1].property;
      }
      else if (j == div2){
	sideinfo.where = EXTERIOR;
	sideinfo.index = exteriorside[2*div1 + div2 - i].index;
	sideinfo.property = exteriorside[2*div1 + div2 - i].property;
      }
      else{
	sideinfo.where = INTERIOR;
	sideinfo.index = (j-1) * div1 + i;
	sideinfo.property = INTERIOR;
      }
    }

    return sideinfo;
}


NODEINFO rectblk::nodeenquiry(unsigned int i, unsigned int j){

    NODEINFO nodeinfo;
    dPOINT p1, p2, p3, p4;
    unsigned long deno;


    if(j == 0){
      nodeinfo.where = EXTERIOR;
      nodeinfo.index = exteriornode[i].index;
      nodeinfo.property = exteriornode[i].property;
      nodeinfo.x = exteriornode[i].x;
      nodeinfo.y = exteriornode[i].y;
    }
    else if (j == div2){
      nodeinfo.where = EXTERIOR;
      nodeinfo.index = exteriornode[2*div1 + div2 - i].index;
      nodeinfo.property = exteriornode[2*div1 + div2 - i].property;
      nodeinfo.x = exteriornode[2*div1 + div2 - i].x;
      nodeinfo.y = exteriornode[2*div1 + div2 - i].y;
    }
    else if (i == 0){
      nodeinfo.where = EXTERIOR;
      nodeinfo.index = exteriornode[2 * (div1 + div2) - j].index;
      nodeinfo.property = exteriornode[2 * (div1 + div2) - j].property;
      nodeinfo.x = exteriornode[2 * (div1 + div2) - j].x;
      nodeinfo.y = exteriornode[2 * (div1 + div2) - j].y;
    }
    else if (i == div1){
      nodeinfo.where = EXTERIOR;
      nodeinfo.index = exteriornode[div1 + j].index;
      nodeinfo.property = exteriornode[div1 + j].property;
      nodeinfo.x = exteriornode[div1 + j].x;
      nodeinfo.y = exteriornode[div1 + j].y;
    }
    else{
      nodeinfo.where = INTERIOR;
      nodeinfo.index = (j - 1) * (div1 - 1) + i;
      nodeinfo.property = INTERIOR;


      vertex(&p1, 1);
      vertex(&p2, 2);
      vertex(&p3, 3);
      vertex(&p4, 4);

      deno = long(div1) * long(div2);

      nodeinfo.x = ((p3.x * (div1 - i) + p4.x * i) * (div2 - j) +
	    (p2.x * (div1 - i) + p1.x * i) * j) / deno; 
      nodeinfo.y = ((p3.y * (div1 - i) + p4.y * i) * (div2 - j) +
	    (p2.y * (div1 - i) + p1.y * i) * j) / deno; 

    }

    return nodeinfo;
}



int quadblk::savegdfinfo(FILE *stream){

   unsigned int i, j,  k, n;

   n = exteriorelmtcount();

   if(fprintf(stream, "%u %u\n", QUADCODE, id) < 0) return FALSE;
   if(fprintf(stream, "%lg %lg %lg %lg %lg %lg\n",
     epsilon, sigma, mu, sigma_mu, xloc, yloc) < 0) return FALSE;
   if(fprintf(stream, "%lg %lg %lg %lg %lg %lg %u %u\n",
     cornerx1, cornery1, cornerx2, cornery2, cornerx3, cornery3, div1, div2)
     < 0) return FALSE;
   for(i = 0; i < n; i++){
     if(fprintf(stream, "%u %u %u %u %u %u %u\n",
       i, exteriornode[i].node, exteriornode[i].cellid, exteriornode[i].property & PROPERTYMASK,
	  exteriorside[i].side, exteriorside[i].cellid, exteriorside[i].property & PROPERTYMASK)
       < 0) return FALSE;
   }
   if(fprintf(stream, "\n") < 0) return FALSE;
   return TRUE;
}

int quadblk::readgdfinfo(FILE *stream){

   unsigned int i, j,  k, n;

   fscanf(stream, "%lf %lf %lf %lf %lf %lf\n",
     &epsilon, &sigma, &mu, &sigma_mu, &xloc, &yloc);
   fscanf(stream, "%lf %lf %lf %lf %lf %lf %u %u\n",
     &cornerx1, &cornery1, &cornerx2, &cornery2, &cornerx3, &cornery3, &div1, &div2);

   n = exteriorelmtcount();

   exteriornode = new surnode[n];
   exteriorside = new surside[n];

   generatenodes();
   for(i = 0; i < n; i++){
     fscanf(stream, "%*u %u %u %u %u %u %u\n",
       &exteriornode[i].node, &exteriornode[i].cellid, &exteriornode[i].property,
       &exteriorside[i].side, &exteriorside[i].cellid, &exteriorside[i].property);
   }
   return TRUE;
}


int triblk::savegdfinfo(FILE *stream){

   unsigned int i, j,  k, n;

   n = exteriorelmtcount();

   if(fprintf(stream, "%u %u\n", TRICODE, id) < 0) return FALSE;
   if(fprintf(stream, "%lg %lg %lg %lg %lg %lg\n",
     epsilon, sigma, mu, sigma_mu, xloc, yloc) < 0) return FALSE;
   if(fprintf(stream, "%lg %lg %lg %lg %u\n",
     cornerx1, cornery1, cornerx2, cornery2, div1) < 0) return FALSE;

   for(i = 0; i < n; i++){
     if(fprintf(stream, "%u %u %u %u %u %u %u\n",
       i, exteriornode[i].node, exteriornode[i].cellid, exteriornode[i].property & PROPERTYMASK,
	  exteriorside[i].side, exteriorside[i].cellid, exteriorside[i].property & PROPERTYMASK)
       < 0) return FALSE;
   }
   if(fprintf(stream, "\n") < 0) return FALSE;
   return TRUE;
}


int triblk::readgdfinfo(FILE *stream){

   unsigned int i, j,  k, n;


   fscanf(stream, "%lf %lf %lf %lf %lf %lf\n",
     &epsilon, &sigma, &mu, &sigma_mu, &xloc, &yloc);
   fscanf(stream, "%lf %lf %lf %lf %u\n",
     &cornerx1, &cornery1, &cornerx2, &cornery2, &div1);

   n = exteriorelmtcount();
   exteriornode = new surnode[n];
   exteriorside = new surside[n];

   generatenodes();

   for(i = 0; i < n; i++){
     fscanf(stream, "%*u %u %u %u %u %u %u\n",
       &exteriornode[i].node, &exteriornode[i].cellid, &exteriornode[i].property,
       &exteriorside[i].side, &exteriorside[i].cellid, &exteriorside[i].property);
   }
   return TRUE;
}


/******************************************************************/
/*                                                                */
/* Member functions for quadblk                                   */
/*                                                                */
/******************************************************************/

// constructor: create objects

quadblk::quadblk(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){
    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;
    ordervertex();
    div1     = d1;
    div2     = d2;
    exteriornode = new surnode[exteriorelmtcount()];
    exteriorside = new surside[exteriorelmtcount()];
    assignid();
    generatenodes();
    initnodes();
    initsides();
}

//  Node assignment:

//   v2  12  11 10  9  8  v1
//       13            7
//       14            6
//       15            5
//   v3   0  1  2   3  4  v4

void quadblk::generatenodes(void){

    int i;
    double sintheta, costheta;
    double x1, y1, x2, y2;
    dPOINT p, p1, p2, p3, p4;

    p.x = cornerx2;
    p.y = cornery2;
    p1.x = xloc + p.x;
    p1.y = yloc + p.y;

    p.x = cornerx3;
    p.y = cornery3;
    p2.x = xloc + p.x;
    p2.y = yloc + p.y;

    p.x = 0.0;
    p.y = 0.0;
    p3.x = xloc + p.x;
    p3.y = yloc + p.y;

    p.x = cornerx1;
    p.y = cornery1;
    p4.x = xloc + p.x;
    p4.y = yloc + p.y;


    for(i = 0; i < div1; i++){
      p.x = (1.0 - (double)i / div1) * p3.x + ((double)i / div1) * p4.x;
      p.y = (1.0 - (double)i / div1) * p3.y + ((double)i / div1) * p4.y;
      exteriornode[i].x = p.x;
      exteriornode[i].y = p.y;
    }

    for(i = 0; i < div2; i++){
      p.x = (1.0 - (double)i / div2) * p4.x + ((double)i / div2) * p1.x;
      p.y = (1.0 - (double)i / div2) * p4.y + ((double)i / div2) * p1.y;
      exteriornode[div1 + i].x = p.x;
      exteriornode[div1 + i].y = p.y;
    }

    for(i = 0; i < div1; i++){
      p.x = (1.0 - (double)i / div1) * p1.x + ((double)i / div1) * p2.x;
      p.y = (1.0 - (double)i / div1) * p1.y + ((double)i / div1) * p2.y;
      exteriornode[div1 + div2 + i].x = p.x;
      exteriornode[div1 + div2 + i].y = p.y;
    }

    for(i = 0; i < div2; i++){