MFEM v2.0
geom.cpp
Go to the documentation of this file.
00001 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
00002 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
00003 // reserved. See file COPYRIGHT for details.
00004 //
00005 // This file is part of the MFEM library. For more information and source code
00006 // availability see http://mfem.googlecode.com.
00007 //
00008 // MFEM is free software; you can redistribute it and/or modify it under the
00009 // terms of the GNU Lesser General Public License (as published by the Free
00010 // Software Foundation) version 2.1 dated February 1999.
00011 
00012 #include "fem.hpp"
00013 
00014 const char *Geometry::Name[NumGeom] =
00015 { "Point", "Segment", "Triangle", "Square", "Tetrahedron", "Cube" };
00016 
00017 Geometry::Geometry()
00018 {
00019    // Vertices for Geometry::POINT
00020    GeomVert[0] = NULL; // No vertices, dimension is 0
00021 
00022    // Vertices for Geometry::SEGMENT
00023    GeomVert[1] = new IntegrationRule(2);
00024 
00025    GeomVert[1]->IntPoint(0).x = 0.0;
00026    GeomVert[1]->IntPoint(1).x = 1.0;
00027 
00028    // Vertices for Geometry::TRIANGLE
00029    GeomVert[2] = new IntegrationRule(3);
00030 
00031    GeomVert[2]->IntPoint(0).x = 0.0;
00032    GeomVert[2]->IntPoint(0).y = 0.0;
00033 
00034    GeomVert[2]->IntPoint(1).x = 1.0;
00035    GeomVert[2]->IntPoint(1).y = 0.0;
00036 
00037    GeomVert[2]->IntPoint(2).x = 0.0;
00038    GeomVert[2]->IntPoint(2).y = 1.0;
00039 
00040    // Vertices for Geometry::SQUARE
00041    GeomVert[3] = new IntegrationRule(4);
00042 
00043    GeomVert[3]->IntPoint(0).x = 0.0;
00044    GeomVert[3]->IntPoint(0).y = 0.0;
00045 
00046    GeomVert[3]->IntPoint(1).x = 1.0;
00047    GeomVert[3]->IntPoint(1).y = 0.0;
00048 
00049    GeomVert[3]->IntPoint(2).x = 1.0;
00050    GeomVert[3]->IntPoint(2).y = 1.0;
00051 
00052    GeomVert[3]->IntPoint(3).x = 0.0;
00053    GeomVert[3]->IntPoint(3).y = 1.0;
00054 
00055    // Vertices for Geometry::TETRAHEDRON
00056    GeomVert[4] = new IntegrationRule(4);
00057    GeomVert[4]->IntPoint(0).x = 0.0;
00058    GeomVert[4]->IntPoint(0).y = 0.0;
00059    GeomVert[4]->IntPoint(0).z = 0.0;
00060 
00061    GeomVert[4]->IntPoint(1).x = 1.0;
00062    GeomVert[4]->IntPoint(1).y = 0.0;
00063    GeomVert[4]->IntPoint(1).z = 0.0;
00064 
00065    GeomVert[4]->IntPoint(2).x = 0.0;
00066    GeomVert[4]->IntPoint(2).y = 1.0;
00067    GeomVert[4]->IntPoint(2).z = 0.0;
00068 
00069    GeomVert[4]->IntPoint(3).x = 0.0;
00070    GeomVert[4]->IntPoint(3).y = 0.0;
00071    GeomVert[4]->IntPoint(3).z = 1.0;
00072 
00073    // Vertices for Geometry::CUBE
00074    GeomVert[5] = new IntegrationRule(8);
00075 
00076    GeomVert[5]->IntPoint(0).x = 0.0;
00077    GeomVert[5]->IntPoint(0).y = 0.0;
00078    GeomVert[5]->IntPoint(0).z = 0.0;
00079 
00080    GeomVert[5]->IntPoint(1).x = 1.0;
00081    GeomVert[5]->IntPoint(1).y = 0.0;
00082    GeomVert[5]->IntPoint(1).z = 0.0;
00083 
00084    GeomVert[5]->IntPoint(2).x = 1.0;
00085    GeomVert[5]->IntPoint(2).y = 1.0;
00086    GeomVert[5]->IntPoint(2).z = 0.0;
00087 
00088    GeomVert[5]->IntPoint(3).x = 0.0;
00089    GeomVert[5]->IntPoint(3).y = 1.0;
00090    GeomVert[5]->IntPoint(3).z = 0.0;
00091 
00092    GeomVert[5]->IntPoint(4).x = 0.0;
00093    GeomVert[5]->IntPoint(4).y = 0.0;
00094    GeomVert[5]->IntPoint(4).z = 1.0;
00095 
00096    GeomVert[5]->IntPoint(5).x = 1.0;
00097    GeomVert[5]->IntPoint(5).y = 0.0;
00098    GeomVert[5]->IntPoint(5).z = 1.0;
00099 
00100    GeomVert[5]->IntPoint(6).x = 1.0;
00101    GeomVert[5]->IntPoint(6).y = 1.0;
00102    GeomVert[5]->IntPoint(6).z = 1.0;
00103 
00104    GeomVert[5]->IntPoint(7).x = 0.0;
00105    GeomVert[5]->IntPoint(7).y = 1.0;
00106    GeomVert[5]->IntPoint(7).z = 1.0;
00107 
00108    GeomCenter[POINT].x = 0.0;
00109    GeomCenter[POINT].y = 0.0;
00110    GeomCenter[POINT].z = 0.0;
00111 
00112    GeomCenter[SEGMENT].x = 0.5;
00113    GeomCenter[SEGMENT].y = 0.0;
00114    GeomCenter[SEGMENT].z = 0.0;
00115 
00116    GeomCenter[TRIANGLE].x = 1.0 / 3.0;
00117    GeomCenter[TRIANGLE].y = 1.0 / 3.0;
00118    GeomCenter[TRIANGLE].z = 0.0;
00119 
00120    GeomCenter[SQUARE].x = 0.5;
00121    GeomCenter[SQUARE].y = 0.5;
00122    GeomCenter[SQUARE].z = 0.0;
00123 
00124    GeomCenter[TETRAHEDRON].x = 0.25;
00125    GeomCenter[TETRAHEDRON].y = 0.25;
00126    GeomCenter[TETRAHEDRON].z = 0.25;
00127 
00128    GeomCenter[CUBE].x = 0.5;
00129    GeomCenter[CUBE].y = 0.5;
00130    GeomCenter[CUBE].z = 0.5;
00131 
00132    PerfGeomToGeomJac[POINT]       = NULL;
00133    PerfGeomToGeomJac[SEGMENT]     = NULL;
00134    PerfGeomToGeomJac[TRIANGLE]    = new DenseMatrix(2);
00135    PerfGeomToGeomJac[SQUARE]      = NULL;
00136    PerfGeomToGeomJac[TETRAHEDRON] = new DenseMatrix(3);
00137    PerfGeomToGeomJac[CUBE]        = NULL;
00138 
00139    {
00140       Linear2DFiniteElement TriFE;
00141       IsoparametricTransformation tri_T;
00142       tri_T.SetFE(&TriFE);
00143       GetPerfPointMat (TRIANGLE, tri_T.GetPointMat());
00144       tri_T.SetIntPoint(&GeomCenter[TRIANGLE]);
00145       CalcInverse(tri_T.Jacobian(), *PerfGeomToGeomJac[TRIANGLE]);
00146    }
00147    {
00148       Linear3DFiniteElement TetFE;
00149       IsoparametricTransformation tet_T;
00150       tet_T.SetFE(&TetFE);
00151       GetPerfPointMat (TETRAHEDRON, tet_T.GetPointMat());
00152       tet_T.SetIntPoint(&GeomCenter[TETRAHEDRON]);
00153       CalcInverse(tet_T.Jacobian(), *PerfGeomToGeomJac[TETRAHEDRON]);
00154    }
00155 }
00156 
00157 Geometry::~Geometry()
00158 {
00159    for (int i = 0; i < NumGeom; i++)
00160    {
00161       delete PerfGeomToGeomJac[i];
00162       delete GeomVert[i];
00163    }
00164 }
00165 
00166 const IntegrationRule * Geometry::GetVertices(int GeomType)
00167 {
00168    switch (GeomType)
00169    {
00170    case Geometry::POINT:       return GeomVert[0];
00171    case Geometry::SEGMENT:     return GeomVert[1];
00172    case Geometry::TRIANGLE:    return GeomVert[2];
00173    case Geometry::SQUARE:      return GeomVert[3];
00174    case Geometry::TETRAHEDRON: return GeomVert[4];
00175    case Geometry::CUBE:        return GeomVert[5];
00176    default:
00177       mfem_error ("Geometry::GetVertices(...)");
00178    }
00179    // make some compilers happy.
00180    return GeomVert[0];
00181 }
00182 
00183 void Geometry::GetPerfPointMat (int GeomType, DenseMatrix &pm)
00184 {
00185    switch (GeomType)
00186    {
00187    case Geometry::SEGMENT:
00188    {
00189       pm.SetSize (1, 2);
00190       pm(0,0) = 0.0;
00191       pm(0,1) = 1.0;
00192    }
00193    break;
00194 
00195    case Geometry::TRIANGLE:
00196    {
00197       pm.SetSize (2, 3);
00198       pm(0,0) = 0.0;  pm(1,0) = 0.0;
00199       pm(0,1) = 1.0;  pm(1,1) = 0.0;
00200       pm(0,2) = 0.5;  pm(1,2) = 0.86602540378443864676;
00201    }
00202    break;
00203 
00204    case Geometry::SQUARE:
00205    {
00206       pm.SetSize (2, 4);
00207       pm(0,0) = 0.0;  pm(1,0) = 0.0;
00208       pm(0,1) = 1.0;  pm(1,1) = 0.0;
00209       pm(0,2) = 1.0;  pm(1,2) = 1.0;
00210       pm(0,3) = 0.0;  pm(1,3) = 1.0;
00211    }
00212    break;
00213 
00214    case Geometry::TETRAHEDRON:
00215    {
00216       pm.SetSize (3, 4);
00217       pm(0,0) = 0.0;  pm(1,0) = 0.0;  pm(2,0) = 0.0;
00218       pm(0,1) = 1.0;  pm(1,1) = 0.0;  pm(2,1) = 0.0;
00219       pm(0,2) = 0.5;  pm(1,2) = 0.86602540378443864676;  pm(2,2) = 0.0;
00220       pm(0,3) = 0.5;  pm(1,3) = 0.28867513459481288225;
00221       pm(2,3) = 0.81649658092772603273;
00222    }
00223    break;
00224 
00225    case Geometry::CUBE:
00226    {
00227       pm.SetSize (3, 8);
00228       pm(0,0) = 0.0;  pm(1,0) = 0.0;  pm(2,0) = 0.0;
00229       pm(0,1) = 1.0;  pm(1,1) = 0.0;  pm(2,1) = 0.0;
00230       pm(0,2) = 1.0;  pm(1,2) = 1.0;  pm(2,2) = 0.0;
00231       pm(0,3) = 0.0;  pm(1,3) = 1.0;  pm(2,3) = 0.0;
00232       pm(0,4) = 0.0;  pm(1,4) = 0.0;  pm(2,4) = 1.0;
00233       pm(0,5) = 1.0;  pm(1,5) = 0.0;  pm(2,5) = 1.0;
00234       pm(0,6) = 1.0;  pm(1,6) = 1.0;  pm(2,6) = 1.0;
00235       pm(0,7) = 0.0;  pm(1,7) = 1.0;  pm(2,7) = 1.0;
00236    }
00237    break;
00238 
00239    default:
00240       mfem_error ("Geometry::GetPerfPointMat (...)");
00241    }
00242 }
00243 
00244 void Geometry::JacToPerfJac(int GeomType, const DenseMatrix &J,
00245                             DenseMatrix &PJ) const
00246 {
00247    if (PerfGeomToGeomJac[GeomType])
00248    {
00249       Mult(J, *PerfGeomToGeomJac[GeomType], PJ);
00250    }
00251    else
00252    {
00253       PJ = J;
00254    }
00255 }
00256 
00257 const int Geometry::NumGeom;
00258 
00259 const int Geometry::NumBdrArray[] = { 0, 2, 3, 4, 4, 6 };
00260 
00261 Geometry Geometries;
00262 
00263 
00264 GeometryRefiner::GeometryRefiner()
00265 {
00266    type = 0;
00267    for (int i = 0; i < Geometry::NumGeom; i++)
00268    {
00269       RGeom[i] = NULL;
00270       IntPts[i] = NULL;
00271    }
00272 }
00273 
00274 GeometryRefiner::~GeometryRefiner()
00275 {
00276    for (int i = 0; i < Geometry::NumGeom; i++)
00277    {
00278       delete RGeom[i];
00279       delete IntPts[i];
00280    }
00281 }
00282 
00283 RefinedGeometry * GeometryRefiner::Refine (int Geom, int Times, int ETimes)
00284 {
00285    int i, j, k, l;
00286 
00287    const double *cp = NULL;
00288    if (type)
00289       cp = poly1d.ClosedPoints(Times);
00290 
00291    switch (Geom)
00292    {
00293    case Geometry::SEGMENT:
00294    {
00295       const int g = Geometry::SEGMENT;
00296       if (RGeom[g] != NULL && RGeom[g]->Times == Times)
00297          return RGeom[g];
00298       delete RGeom[g];
00299       RGeom[g] = new RefinedGeometry(Times+1, 2*Times, 0);
00300       RGeom[g]->Times = Times;
00301       RGeom[g]->ETimes = 0;
00302       for (i = 0; i <= Times; i++)
00303       {
00304          IntegrationPoint &ip = RGeom[g]->RefPts.IntPoint(i);
00305          ip.x = (type == 0) ? double(i) / Times : cp[i];
00306       }
00307       Array<int> &G = RGeom[g]->RefGeoms;
00308       for (i = 0; i < Times; i++)
00309       {
00310          G[2*i+0] = i;
00311          G[2*i+1] = i+1;
00312       }
00313 
00314       return RGeom[g];
00315    }
00316 
00317    case Geometry::TRIANGLE:
00318    {
00319       if (RGeom[2] != NULL && RGeom[2]->Times == Times &&
00320           RGeom[2]->ETimes == ETimes)
00321          return RGeom[2];
00322 
00323       if (RGeom[2] != NULL)
00324          delete RGeom[2];
00325       RGeom[2] = new RefinedGeometry ((Times+1)*(Times+2)/2, 3*Times*Times,
00326                                       3*Times*(ETimes+1));
00327       RGeom[2]->Times = Times;
00328       RGeom[2]->ETimes = ETimes;
00329       for (k = j = 0; j <= Times; j++)
00330          for (i = 0; i <= Times-j; i++, k++)
00331          {
00332             IntegrationPoint &ip = RGeom[2]->RefPts.IntPoint(k);
00333             if (type == 0)
00334             {
00335                ip.x = double(i) / Times;
00336                ip.y = double(j) / Times;
00337             }
00338             else
00339             {
00340                ip.x = cp[i]/(cp[i] + cp[j] + cp[Times-i-j]);
00341                ip.y = cp[j]/(cp[i] + cp[j] + cp[Times-i-j]);
00342             }
00343          }
00344       Array<int> &G = RGeom[2]->RefGeoms;
00345       for (l = k = j = 0; j < Times; j++, k++)
00346          for (i = 0; i < Times-j; i++, k++)
00347          {
00348             G[l++] = k;
00349             G[l++] = k+1;
00350             G[l++] = k+Times-j+1;
00351             if (i+j+1 < Times)
00352             {
00353                G[l++] = k+1;
00354                G[l++] = k+Times-j+2;
00355                G[l++] = k+Times-j+1;
00356             }
00357          }
00358       Array<int> &E = RGeom[2]->RefEdges;
00359       // horizontal edges
00360       for (l = k = 0; k < Times; k += Times/ETimes)
00361       {
00362          j = k*(Times+1)-((k-1)*k)/2;
00363          for (i = 0; i < Times-k; i++)
00364          {
00365             E[l++] = j; j++;
00366             E[l++] = j;
00367          }
00368       }
00369       // diagonal edges
00370       for (k = Times; k > 0; k -= Times/ETimes)
00371       {
00372          j = k;
00373          for (i = 0; i < k; i++)
00374          {
00375             E[l++] = j; j += Times-i;
00376             E[l++] = j;
00377          }
00378       }
00379       // vertical edges
00380       for (k = 0; k < Times; k += Times/ETimes)
00381       {
00382          j = k;
00383          for (i = 0; i < Times-k; i++)
00384          {
00385             E[l++] = j; j += Times-i+1;
00386             E[l++] = j;
00387          }
00388       }
00389 
00390       return RGeom[2];
00391    }
00392 
00393    case Geometry::SQUARE:
00394    {
00395       if (RGeom[3] != NULL && RGeom[3]->Times == Times &&
00396           RGeom[3]->ETimes == ETimes)
00397          return RGeom[3];
00398 
00399       if (RGeom[3] != NULL)
00400          delete RGeom[3];
00401       RGeom[3] = new RefinedGeometry ((Times+1)*(Times+1), 4*Times*Times,
00402                                       4*(ETimes+1)*Times);
00403       RGeom[3]->Times = Times;
00404       RGeom[3]->ETimes = ETimes;
00405       for (k = j = 0; j <= Times; j++)
00406          for (i = 0; i <= Times; i++, k++)
00407          {
00408             IntegrationPoint &ip = RGeom[3]->RefPts.IntPoint(k);
00409             if (type == 0)
00410             {
00411                ip.x = double(i) / Times;
00412                ip.y = double(j) / Times;
00413             }
00414             else
00415             {
00416                ip.x = cp[i];
00417                ip.y = cp[j];
00418             }
00419          }
00420       Array<int> &G = RGeom[3]->RefGeoms;
00421       for (l = k = j = 0; j < Times; j++, k++)
00422          for (i = 0; i < Times; i++, k++)
00423          {
00424             G[l++] = k;
00425             G[l++] = k+1;
00426             G[l++] = k+Times+2;
00427             G[l++] = k+Times+1;
00428          }
00429       Array<int> &E = RGeom[3]->RefEdges;
00430       // horizontal edges
00431       for (l = k = 0; k <= Times; k += Times/ETimes)
00432       {
00433          for (i = 0, j = k*(Times+1); i < Times; i++)
00434          {
00435             E[l++] = j; j++;
00436             E[l++] = j;
00437          }
00438       }
00439       // vertical edges (in right-to-left order)
00440       for (k = Times; k >= 0; k -= Times/ETimes)
00441       {
00442          for (i = 0, j = k; i < Times; i++)
00443          {
00444             E[l++] = j; j += Times+1;
00445             E[l++] = j;
00446          }
00447       }
00448 
00449       return RGeom[3];
00450    }
00451 
00452    case Geometry::CUBE:
00453    {
00454       const int g = Geometry::CUBE;
00455       if (RGeom[g] != NULL && RGeom[g]->Times == Times &&
00456           RGeom[g]->ETimes == ETimes)
00457          return RGeom[g];
00458 
00459       if (RGeom[g] != NULL)
00460          delete RGeom[g];
00461       RGeom[g] = new RefinedGeometry ((Times+1)*(Times+1)*(Times+1),
00462                                       8*Times*Times*Times, 0);
00463       RGeom[g]->Times = Times;
00464       RGeom[g]->ETimes = ETimes;
00465       for (l = k = 0; k <= Times; k++)
00466          for (j = 0; j <= Times; j++)
00467             for (i = 0; i <= Times; i++, l++)
00468             {
00469                IntegrationPoint &ip = RGeom[g]->RefPts.IntPoint(l);
00470                if (type == 0)
00471                {
00472                   ip.x = double(i) / Times;
00473                   ip.y = double(j) / Times;
00474                   ip.z = double(k) / Times;
00475                }
00476                else
00477                {
00478                   ip.x = cp[i];
00479                   ip.y = cp[j];
00480                   ip.z = cp[k];
00481                }
00482             }
00483       Array<int> &G = RGeom[g]->RefGeoms;
00484       for (l = k = 0; k < Times; k++)
00485          for (j = 0; j < Times; j++)
00486             for (i = 0; i < Times; i++)
00487             {
00488                G[l++] = i+0 + (j+0 + (k+0) * (Times+1)) * (Times+1);
00489                G[l++] = i+1 + (j+0 + (k+0) * (Times+1)) * (Times+1);
00490                G[l++] = i+1 + (j+1 + (k+0) * (Times+1)) * (Times+1);
00491                G[l++] = i+0 + (j+1 + (k+0) * (Times+1)) * (Times+1);
00492                G[l++] = i+0 + (j+0 + (k+1) * (Times+1)) * (Times+1);
00493                G[l++] = i+1 + (j+0 + (k+1) * (Times+1)) * (Times+1);
00494                G[l++] = i+1 + (j+1 + (k+1) * (Times+1)) * (Times+1);
00495                G[l++] = i+0 + (j+1 + (k+1) * (Times+1)) * (Times+1);
00496             }
00497 
00498       return RGeom[g];
00499    }
00500 
00501    case Geometry::TETRAHEDRON:
00502    {
00503       const int g = Geometry::TETRAHEDRON;
00504       if (RGeom[g] != NULL && RGeom[g]->Times == Times &&
00505           RGeom[g]->ETimes == ETimes)
00506          return RGeom[g];
00507 
00508       if (RGeom[g] != NULL)
00509          delete RGeom[g];
00510 
00511       // subdivide the tetrahedron with vertices
00512       // (0,0,0), (0,0,1), (1,1,1), (0,1,1)
00513 
00514       // vertices: 0 <= i <= j <= k <= Times
00515       // (3-combination with repetitions)
00516       // number of vertices: (n+3)*(n+2)*(n+1)/6, n = Times
00517 
00518       // elements: the vertices are: v1=(i,j,k), v2=v1+u1, v3=v2+u2, v4=v3+u3
00519       // where 0 <= i <= j <= k <= n-1 and
00520       // u1,u2,u3 is a permutation of (1,0,0),(0,1,0),(0,0,1)
00521       // such that all v2,v3,v4 have non-decreasing components
00522       // number of elements: n^3
00523 
00524       const int n = Times;
00525       RGeom[g] = new RefinedGeometry((n+3)*(n+2)*(n+1)/6, 4*n*n*n, 0);
00526       // enumerate and define the vertices
00527       Array<int> vi((n+1)*(n+1)*(n+1));
00528       vi = -1;
00529       int m = 0;
00530       for (k = 0; k <= n; k++)
00531          for (j = 0; j <= k; j++)
00532             for (i = 0; i <= j; i++)
00533             {
00534                IntegrationPoint &ip = RGeom[g]->RefPts.IntPoint(m);
00535                // map the coordinates to the reference tetrahedron
00536                // (0,0,0) -> (0,0,0)
00537                // (0,0,1) -> (1,0,0)
00538                // (1,1,1) -> (0,1,0)
00539                // (0,1,1) -> (0,0,1)
00540                if (type == 0)
00541                {
00542                   ip.x = double(k - j) / n;
00543                   ip.y = double(i) / n;
00544                   ip.z = double(j - i) / n;
00545                }
00546                else
00547                {
00548                   double w = cp[k-j] + cp[i] + cp[j-i] + cp[Times-k];
00549                   ip.x = cp[k-j]/w;
00550                   ip.y = cp[i]/w;
00551                   ip.z = cp[j-i]/w;
00552                }
00553                l = i + (j + k * (n+1)) * (n+1);
00554                vi[l] = m;
00555                m++;
00556             }
00557       if (m != (n+3)*(n+2)*(n+1)/6)
00558          mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #1");
00559       // elements
00560       Array<int> &G = RGeom[g]->RefGeoms;
00561       m = 0;
00562       for (k = 0; k < n; k++)
00563          for (j = 0; j <= k; j++)
00564             for (i = 0; i <= j; i++)
00565             {
00566                // the ordering of the vertices is chosen to ensure:
00567                // 1) correct orientation
00568                // 2) the x,y,z edges are in the set of edges
00569                //    {(0,1),(2,3), (0,2),(1,3)}
00570                //    (goal is to ensure that subsequent refinement using
00571                //    this procedure preserves the six tetrahedral shapes)
00572 
00573                // zyx: (i,j,k)-(i,j,k+1)-(i+1,j+1,k+1)-(i,j+1,k+1)
00574                G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
00575                G[m++] = vi[i+0 + (j+0 + (k+1) * (n+1)) * (n+1)];
00576                G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
00577                G[m++] = vi[i+0 + (j+1 + (k+1) * (n+1)) * (n+1)];
00578                if (j < k)
00579                {
00580                   // yzx: (i,j,k)-(i+1,j+1,k+1)-(i,j+1,k)-(i,j+1,k+1)
00581                   G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
00582                   G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
00583                   G[m++] = vi[i+0 + (j+1 + (k+0) * (n+1)) * (n+1)];
00584                   G[m++] = vi[i+0 + (j+1 + (k+1) * (n+1)) * (n+1)];
00585                   // yxz: (i,j,k)-(i,j+1,k)-(i+1,j+1,k+1)-(i+1,j+1,k)
00586                   G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
00587                   G[m++] = vi[i+0 + (j+1 + (k+0) * (n+1)) * (n+1)];
00588                   G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
00589                   G[m++] = vi[i+1 + (j+1 + (k+0) * (n+1)) * (n+1)];
00590                }
00591                if (i < j)
00592                {
00593                   // xzy: (i,j,k)-(i+1,j,k)-(i+1,j+1,k+1)-(i+1,j,k+1)
00594                   G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
00595                   G[m++] = vi[i+1 + (j+0 + (k+0) * (n+1)) * (n+1)];
00596                   G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
00597                   G[m++] = vi[i+1 + (j+0 + (k+1) * (n+1)) * (n+1)];
00598                   if (j < k)
00599                   {
00600                      // xyz: (i,j,k)-(i+1,j+1,k+1)-(i+1,j,k)-(i+1,j+1,k)
00601                      G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
00602                      G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
00603                      G[m++] = vi[i+1 + (j+0 + (k+0) * (n+1)) * (n+1)];
00604                      G[m++] = vi[i+1 + (j+1 + (k+0) * (n+1)) * (n+1)];
00605                   }
00606                   // zxy: (i,j,k)-(i+1,j+1,k+1)-(i,j,k+1)-(i+1,j,k+1)
00607                   G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
00608                   G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
00609                   G[m++] = vi[i+0 + (j+0 + (k+1) * (n+1)) * (n+1)];
00610                   G[m++] = vi[i+1 + (j+0 + (k+1) * (n+1)) * (n+1)];
00611                }
00612             }
00613       if (m != 4*n*n*n)
00614          mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #2");
00615       for (i = 0; i < m; i++)
00616          if (G[i] < 0)
00617             mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #3");
00618 
00619       return RGeom[g];
00620    }
00621 
00622    default:
00623 
00624       return RGeom[0];
00625    }
00626 }
00627 
00628 const IntegrationRule *GeometryRefiner::RefineInterior(int Geom, int Times)
00629 {
00630    int g = Geom;
00631 
00632    switch (g)
00633    {
00634    case Geometry::SEGMENT:
00635    {
00636       if (Times < 2)
00637          return NULL;
00638       if (IntPts[g] == NULL || IntPts[g]->GetNPoints() != Times-1)
00639       {
00640          delete IntPts[g];
00641          IntPts[g] = new IntegrationRule(Times-1);
00642          for (int i = 1; i < Times; i++)
00643          {
00644             IntegrationPoint &ip = IntPts[g]->IntPoint(i-1);
00645             ip.x = double(i) / Times;
00646             ip.y = ip.z = 0.0;
00647          }
00648       }
00649    }
00650    break;
00651 
00652    case Geometry::TRIANGLE:
00653    {
00654       if (Times < 3)
00655          return NULL;
00656       if (IntPts[g] == NULL ||
00657           IntPts[g]->GetNPoints() != ((Times-1)*(Times-2))/2)
00658       {
00659          delete IntPts[g];
00660          IntPts[g] = new IntegrationRule(((Times-1)*(Times-2))/2);
00661          for (int k = 0, j = 1; j < Times-1; j++)
00662             for (int i = 1; i < Times-j; i++, k++)
00663             {
00664                IntegrationPoint &ip = IntPts[g]->IntPoint(k);
00665                ip.x = double(i) / Times;
00666                ip.y = double(j) / Times;
00667                ip.z = 0.0;
00668             }
00669       }
00670    }
00671    break;
00672 
00673    case Geometry::SQUARE:
00674    {
00675       if (Times < 2)
00676          return NULL;
00677       if (IntPts[g] == NULL || IntPts[g]->GetNPoints() != (Times-1)*(Times-1))
00678       {
00679          delete IntPts[g];
00680          IntPts[g] = new IntegrationRule((Times-1)*(Times-1));
00681          for (int k = 0, j = 1; j < Times; j++)
00682             for (int i = 1; i < Times; i++, k++)
00683             {
00684                IntegrationPoint &ip = IntPts[g]->IntPoint(k);
00685                ip.x = double(i) / Times;
00686                ip.y = double(j) / Times;
00687                ip.z = 0.0;
00688             }
00689       }
00690    }
00691    break;
00692 
00693    default:
00694       mfem_error("GeometryRefiner::RefineInterior(...)");
00695    }
00696 
00697    return IntPts[g];
00698 }
00699 
00700 GeometryRefiner GlobGeometryRefiner;
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines