MFEM v2.0
fespace.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 // Implementation of FiniteElementSpace
00013 
00014 #include "fem.hpp"
00015 #include <math.h>
00016 
00017 int FiniteElementSpace::GetOrder(int i) const
00018 {
00019    int GeomType = mesh->GetElementBaseGeometry(i);
00020    return fec->FiniteElementForGeometry(GeomType)->GetOrder();
00021 }
00022 
00023 
00024 void FiniteElementSpace::DofsToVDofs (Array<int> &dofs) const
00025 {
00026    int i, j, size;
00027 
00028    if (vdim == 1)  return;
00029 
00030    size = dofs.Size();
00031    dofs.SetSize (size * vdim);
00032 
00033    switch(ordering)
00034    {
00035    case Ordering::byNODES:
00036       for (i = 1; i < vdim; i++)
00037          for (j = 0; j < size; j++)
00038             if (dofs[j] < 0)
00039                dofs[size * i + j] = -1 - ( ndofs * i + (-1-dofs[j]) );
00040             else
00041                dofs[size * i + j] = ndofs * i + dofs[j];
00042       break;
00043 
00044    case Ordering::byVDIM:
00045       for (i = vdim-1; i >= 0; i--)
00046          for (j = 0; j < size; j++)
00047             if (dofs[j] < 0)
00048                dofs[size * i + j] = -1 - ( (-1-dofs[j]) * vdim + i );
00049             else
00050                dofs[size * i + j] = dofs[j] * vdim + i;
00051       break;
00052    }
00053 }
00054 
00055 void FiniteElementSpace::DofsToVDofs(int vd, Array<int> &dofs) const
00056 {
00057    if (vdim == 1)
00058       return;
00059    if (ordering == Ordering::byNODES)
00060    {
00061       for (int i = 0; i < dofs.Size(); i++)
00062       {
00063          int dof = dofs[i];
00064          if (dof < 0)
00065             dofs[i] = -1 - ((-1-dof) + vd * ndofs);
00066          else
00067             dofs[i] = dof + vd * ndofs;
00068       }
00069    }
00070    else
00071    {
00072       for (int i = 0; i < dofs.Size(); i++)
00073       {
00074          int dof = dofs[i];
00075          if (dof < 0)
00076             dofs[i] = -1 - ((-1-dof) * vdim + vd);
00077          else
00078             dofs[i] = dof * vdim + vd;
00079       }
00080    }
00081 }
00082 
00083 int FiniteElementSpace::DofToVDof (int dof, int vd) const
00084 {
00085    if (vdim == 1)
00086       return dof;
00087    if (ordering == Ordering::byNODES)
00088    {
00089       if (dof < 0)
00090          return -1 - ((-1-dof) + vd * ndofs);
00091       else
00092          return dof + vd * ndofs;
00093    }
00094    if (dof < 0)
00095       return -1 - ((-1-dof) * vdim + vd);
00096    else
00097       return dof * vdim + vd;
00098 }
00099 
00100 // static function
00101 void FiniteElementSpace::AdjustVDofs (Array<int> &vdofs)
00102 {
00103    int n = vdofs.Size(), *vdof = vdofs;
00104    for (int i = 0; i < n; i++)
00105    {
00106       int j;
00107       if ((j=vdof[i]) < 0)
00108          vdof[i] = -1-j;
00109    }
00110 }
00111 
00112 void FiniteElementSpace::GetElementVDofs(int iE, Array<int> &dofs) const
00113 {
00114    GetElementDofs(iE, dofs);
00115    DofsToVDofs (dofs);
00116 }
00117 
00118 void FiniteElementSpace::GetBdrElementVDofs (int iE, Array<int> &dofs) const
00119 {
00120    GetBdrElementDofs(iE, dofs);
00121    DofsToVDofs (dofs);
00122 }
00123 
00124 void FiniteElementSpace::GetFaceVDofs (int iF, Array<int> &dofs) const
00125 {
00126    GetFaceDofs (iF, dofs);
00127    DofsToVDofs (dofs);
00128 }
00129 
00130 void FiniteElementSpace::GetEdgeVDofs (int iE, Array<int> &dofs) const
00131 {
00132    GetEdgeDofs (iE, dofs);
00133    DofsToVDofs (dofs);
00134 }
00135 
00136 void FiniteElementSpace::GetElementInteriorVDofs (int i, Array<int> &vdofs)
00137    const
00138 {
00139    GetElementInteriorDofs (i, vdofs);
00140    DofsToVDofs (vdofs);
00141 }
00142 
00143 void FiniteElementSpace::GetEdgeInteriorVDofs (int i, Array<int> &vdofs)
00144    const
00145 {
00146    GetEdgeInteriorDofs (i, vdofs);
00147    DofsToVDofs (vdofs);
00148 }
00149 
00150 void FiniteElementSpace::BuildElementToDofTable()
00151 {
00152    if (elem_dof)
00153       return;
00154 
00155    Table *el_dof = new Table;
00156    Array<int> dofs;
00157    el_dof -> MakeI (mesh -> GetNE());
00158    for (int i = 0; i < mesh -> GetNE(); i++)
00159    {
00160       GetElementDofs (i, dofs);
00161       el_dof -> AddColumnsInRow (i, dofs.Size());
00162    }
00163    el_dof -> MakeJ();
00164    for (int i = 0; i < mesh -> GetNE(); i++)
00165    {
00166       GetElementDofs (i, dofs);
00167       el_dof -> AddConnections (i, (int *)dofs, dofs.Size());
00168    }
00169    el_dof -> ShiftUpI();
00170    elem_dof = el_dof;
00171 }
00172 
00173 void FiniteElementSpace::BuildDofToArrays()
00174 {
00175    if (dof_elem_array.Size())
00176       return;
00177    BuildElementToDofTable();
00178 
00179    dof_elem_array.SetSize (ndofs);
00180    dof_ldof_array.SetSize (ndofs);
00181    dof_elem_array = -1;
00182    for (int i = 0; i < mesh -> GetNE(); i++)
00183    {
00184       const int *dofs = elem_dof -> GetRow(i);
00185       const int n = elem_dof -> RowSize(i);
00186       for (int j = 0; j < n; j++)
00187          if (dof_elem_array[dofs[j]] < 0)
00188          {
00189             dof_elem_array[dofs[j]] = i;
00190             dof_ldof_array[dofs[j]] = j;
00191          }
00192    }
00193 }
00194 
00195 DenseMatrix * FiniteElementSpace::LocalInterpolation
00196 (int k, int num_c_dofs, RefinementType type, Array<int> &rows)
00197 {
00198    int i,j,l=-1;
00199 
00200    // generate refinement data if necessary
00201    for (i = 0; i < RefData.Size(); i++)
00202       if (RefData[i] -> type == type)
00203       { l = i; break; }
00204 
00205    if (l == -1)
00206    { ConstructRefinementData (k,num_c_dofs,type);  l = i; }
00207 
00208    // determine the global indices of the fine dofs on the coarse element
00209    Array<int> l_dofs, g_dofs;
00210    int num_fine_elems = RefData[l] -> num_fine_elems;
00211    rows.SetSize(RefData[l] -> num_fine_dofs);
00212 
00213    for (i = 0; i < num_fine_elems; i++) {
00214       GetElementDofs(mesh->GetFineElem(k,i),g_dofs); // TWO_LEVEL_FINE
00215       RefData[l] -> fl_to_fc -> GetRow (i, l_dofs);
00216       for (j = 0; j < l_dofs.Size(); j++)
00217          rows[l_dofs[j]] = g_dofs[j];
00218    }
00219 
00220    return RefData[l] -> I;
00221 }
00222 
00223 SparseMatrix * FiniteElementSpace::GlobalRestrictionMatrix
00224 (FiniteElementSpace * cfes, int one_vdim)
00225 {
00226    int k;
00227    SparseMatrix * R;
00228    DenseMatrix * r;
00229    Array<int> rows, cols, rs, cs;
00230 
00231    if (one_vdim == -1)
00232       one_vdim = (ordering == Ordering::byNODES) ? 1 : 0;
00233 
00234    mesh -> SetState (Mesh::TWO_LEVEL_COARSE);
00235    if (vdim == 1 || one_vdim == 1)
00236       R = new SparseMatrix (cfes -> GetNDofs(), ndofs);
00237    else
00238       R = new SparseMatrix (cfes -> GetVSize(), vdim*ndofs);
00239 
00240    for (k = 0; k < mesh -> GetNE(); k++)
00241    {
00242       cfes -> GetElementDofs (k, rows);
00243       mesh -> SetState (Mesh::TWO_LEVEL_FINE);
00244       r = LocalInterpolation (k, rows.Size(),
00245                               mesh -> GetRefinementType(k), cols);
00246       if (vdim == 1 || one_vdim == 1)
00247          R -> SetSubMatrixTranspose (rows, cols, *r, 1);
00248       else
00249       {
00250          int d, nr = rows.Size(), nc = cols.Size();
00251          cfes -> DofsToVDofs (rows);
00252          DofsToVDofs (cols);
00253          for (d = 0; d < vdim; d++)
00254          {
00255             rows.GetSubArray (d*nr, nr, rs);
00256             cols.GetSubArray (d*nc, nc, cs);
00257             R -> SetSubMatrixTranspose (rs, cs, *r, 1);
00258          }
00259       }
00260       mesh -> SetState (Mesh::TWO_LEVEL_COARSE);
00261    }
00262 
00263    return R;
00264 }
00265 
00266 void FiniteElementSpace::GetEssentialVDofs (Array<int> &bdr_attr_is_ess,
00267                                             Array<int> &ess_dofs)
00268 {
00269    int i, j, k;
00270    Array<int> vdofs;
00271 
00272    ess_dofs.SetSize (GetVSize());
00273    ess_dofs = 0;
00274 
00275    for (i = 0; i < GetNBE(); i++)
00276       if (bdr_attr_is_ess[GetBdrAttribute (i)-1])
00277       {
00278          GetBdrElementVDofs (i, vdofs);
00279          for (j = 0; j < vdofs.Size(); j++)
00280             if ( (k = vdofs[j]) >= 0 )
00281                ess_dofs[k] = -1;
00282             else
00283                ess_dofs[-1-k] = -1;
00284       }
00285 }
00286 
00287 void FiniteElementSpace::EliminateEssentialBCFromGRM
00288 (FiniteElementSpace *cfes, Array<int> &bdr_attr_is_ess, SparseMatrix *R)
00289 {
00290    int i, j, k, one_vdim;
00291    Array<int> dofs;
00292 
00293    one_vdim = (cfes -> GetNDofs() == R -> Size()) ? 1 : 0;
00294 
00295    mesh -> SetState (Mesh::TWO_LEVEL_COARSE);
00296    if (bdr_attr_is_ess.Size() != 0)
00297       for (i=0; i < cfes -> GetNBE(); i++)
00298          if (bdr_attr_is_ess[cfes -> GetBdrAttribute (i)-1])
00299          {
00300             if (one_vdim == 1)
00301                cfes -> GetBdrElementDofs (i, dofs);
00302             else
00303                cfes -> GetBdrElementVDofs (i, dofs);
00304             for (j=0; j < dofs.Size(); j++)
00305                if ( (k = dofs[j]) >= 0 )
00306                   R -> EliminateRow(k);
00307                else
00308                   R -> EliminateRow(-1-k);
00309          }
00310    R -> Finalize();
00311 }
00312 
00313 SparseMatrix * FiniteElementSpace::GlobalRestrictionMatrix
00314 (FiniteElementSpace * cfes, Array<int> &bdr_attr_is_ess, int one_vdim)
00315 {
00316    SparseMatrix * R;
00317 
00318    R = GlobalRestrictionMatrix (cfes, one_vdim);
00319    EliminateEssentialBCFromGRM (cfes, bdr_attr_is_ess, R);
00320 
00321    return R;
00322 }
00323 
00324 SparseMatrix *
00325 FiniteElementSpace::D2C_GlobalRestrictionMatrix (FiniteElementSpace *cfes)
00326 {
00327    int i, j;
00328    Array<int> d_vdofs, c_vdofs;
00329    SparseMatrix *R;
00330 
00331    R = new SparseMatrix (cfes -> GetVSize(), GetVSize());
00332 
00333    for (i = 0; i < mesh -> GetNE(); i++)
00334    {
00335       this -> GetElementVDofs (i, d_vdofs);
00336       cfes -> GetElementVDofs (i, c_vdofs);
00337 
00338 #ifdef MFEM_DEBUG
00339       if (d_vdofs.Size() != c_vdofs.Size())
00340          mfem_error ("FiniteElementSpace::D2C_GlobalRestrictionMatrix (...)");
00341 #endif
00342 
00343       for (j = 0; j < d_vdofs.Size(); j++)
00344          R -> Set (c_vdofs[j], d_vdofs[j], 1.0);
00345    }
00346 
00347    R -> Finalize();
00348 
00349    return R;
00350 }
00351 
00352 SparseMatrix *
00353 FiniteElementSpace::D2Const_GlobalRestrictionMatrix(FiniteElementSpace *cfes)
00354 {
00355    int i, j;
00356    Array<int> d_dofs, c_dofs;
00357    SparseMatrix *R;
00358 
00359    R = new SparseMatrix (cfes -> GetNDofs(), ndofs);
00360 
00361    for (i = 0; i < mesh -> GetNE(); i++)
00362    {
00363       this -> GetElementDofs (i, d_dofs);
00364       cfes -> GetElementDofs (i, c_dofs);
00365 
00366 #ifdef MFEM_DEBUG
00367       if (c_dofs.Size() != 1)
00368          mfem_error ("FiniteElementSpace::"
00369                      "D2Const_GlobalRestrictionMatrix (...)");
00370 #endif
00371 
00372       for (j = 0; j < d_dofs.Size(); j++)
00373          R -> Set (c_dofs[0], d_dofs[j], 1.0);
00374    }
00375 
00376    R -> Finalize();
00377 
00378    return R;
00379 }
00380 
00381 SparseMatrix *
00382 FiniteElementSpace::H2L_GlobalRestrictionMatrix (FiniteElementSpace *lfes)
00383 {
00384    int i, j;
00385    SparseMatrix *R;
00386    DenseMatrix loc_restr;
00387    Vector shape;
00388    Array<int> l_dofs, h_dofs;
00389 
00390    R = new SparseMatrix (lfes -> GetNDofs(), ndofs);
00391 
00392    const FiniteElement *h_fe = this -> GetFE (0);
00393    const FiniteElement *l_fe = lfes -> GetFE (0);
00394    shape.SetSize (l_fe -> GetDof());
00395    loc_restr.SetSize (l_fe -> GetDof(), h_fe -> GetDof());
00396    for (i = 0; i < h_fe -> GetDof(); i++)
00397    {
00398       l_fe -> CalcShape (h_fe -> GetNodes().IntPoint(i), shape);
00399       for (j = 0; j < shape.Size(); j++)
00400       {
00401          if (fabs (shape(j)) < 1.0e-6)
00402             shape(j) = 0.0;
00403          loc_restr(j,i) = shape(j);
00404       }
00405    }
00406 
00407    for (i = 0; i < mesh -> GetNE(); i++)
00408    {
00409       this -> GetElementDofs (i, h_dofs);
00410       lfes -> GetElementDofs (i, l_dofs);
00411 
00412       R -> SetSubMatrix (l_dofs, h_dofs, loc_restr, 1);
00413    }
00414 
00415    R -> Finalize();
00416 
00417    return R;
00418 }
00419 
00420 FiniteElementSpace::FiniteElementSpace(FiniteElementSpace &fes)
00421 {
00422    mesh = fes.mesh;
00423    vdim = fes.vdim;
00424    ndofs = fes.ndofs;
00425    ordering = fes.ordering;
00426    fec = fes.fec;
00427    nvdofs = fes.nvdofs;
00428    nedofs = fes.nedofs;
00429    nfdofs = fes.nfdofs;
00430    nbdofs = fes.nbdofs;
00431    fdofs = fes.fdofs;
00432    bdofs = fes.bdofs;
00433    // keep 'RefData' in 'fes'
00434    elem_dof = fes.elem_dof;
00435    bdrElem_dof = fes.bdrElem_dof;
00436    Swap(dof_elem_array, fes.dof_elem_array);
00437    Swap(dof_ldof_array, fes.dof_ldof_array);
00438 
00439    NURBSext = fes.NURBSext;
00440    own_ext = 0;
00441 
00442    fes.bdofs = NULL;
00443    fes.fdofs = NULL;
00444    fes.elem_dof = NULL;
00445    fes.bdrElem_dof = NULL;
00446 }
00447 
00448 FiniteElementSpace::FiniteElementSpace(Mesh *m,
00449                                        const FiniteElementCollection *f,
00450                                        int dim, int order)
00451 {
00452    mesh = m;
00453    fec = f;
00454    vdim = dim;
00455    ordering = order;
00456 
00457    const NURBSFECollection *nurbs_fec =
00458       dynamic_cast<const NURBSFECollection *>(fec);
00459    if (nurbs_fec)
00460    {
00461       if (!mesh->NURBSext)
00462       {
00463          mfem_error("FiniteElementSpace::FiniteElementSpace :\n"
00464                     "   NURBS FE space requires NURBS mesh.");
00465       }
00466       else
00467       {
00468          int Order = nurbs_fec->GetOrder();
00469          if (mesh->NURBSext->GetOrder() == Order)
00470          {
00471             NURBSext = mesh->NURBSext;
00472             own_ext = 0;
00473          }
00474          else
00475          {
00476             NURBSext = new NURBSExtension(mesh->NURBSext, Order);
00477             own_ext = 1;
00478          }
00479          UpdateNURBS();
00480       }
00481    }
00482    else
00483    {
00484       NURBSext = NULL;
00485       own_ext = 0;
00486       Constructor();
00487    }
00488 }
00489 
00490 void FiniteElementSpace::UpdateNURBS()
00491 {
00492    nvdofs = 0;
00493    nedofs = 0;
00494    nfdofs = 0;
00495    nbdofs = 0;
00496    fdofs = NULL;
00497    bdofs = NULL;
00498 
00499    dynamic_cast<const NURBSFECollection *>(fec)->Reset();
00500 
00501    ndofs = NURBSext->GetNDof();
00502 
00503    elem_dof = NURBSext->GetElementDofTable();
00504 
00505    bdrElem_dof = NURBSext->GetBdrElementDofTable();
00506 }
00507 
00508 void FiniteElementSpace::Constructor()
00509 {
00510    int i;
00511 
00512    elem_dof = NULL;
00513    bdrElem_dof = NULL;
00514 
00515    nvdofs = mesh->GetNV() * fec->DofForGeometry(Geometry::POINT);
00516 
00517    if ( mesh->Dimension() > 1 )
00518       nedofs = mesh->GetNEdges() * fec->DofForGeometry(Geometry::SEGMENT);
00519    else
00520       nedofs = 0;
00521 
00522    nfdofs = 0;
00523    fdofs = NULL;
00524    if (mesh->Dimension() == 3)
00525    {
00526       // Here we assume that all faces in the mesh have the same base
00527       // geometry -- the base geometry of the 0-th face element.
00528       // The class Mesh assumes the same inside GetFaceBaseGeometry(...).
00529       // Thus we do not need to generate all the faces in the mesh
00530       // if we do not need them.
00531       int fdof = fec->DofForGeometry(mesh->GetFaceBaseGeometry(0));
00532       if (fdof > 0)
00533       {
00534          fdofs = new int[mesh->GetNFaces()+1];
00535          fdofs[0] = 0;
00536          for (i = 0; i < mesh->GetNFaces(); i++)
00537          {
00538             nfdofs += fdof;
00539             // nfdofs += fec->DofForGeometry(mesh->GetFaceBaseGeometry(i));
00540             fdofs[i+1] = nfdofs;
00541          }
00542       }
00543    }
00544 
00545    nbdofs = 0;
00546    bdofs = new int[mesh->GetNE()+1];
00547    bdofs[0] = 0;
00548    for (i = 0; i < mesh->GetNE(); i++)
00549    {
00550       nbdofs += fec->DofForGeometry(mesh->GetElementBaseGeometry(i));
00551       bdofs[i+1] = nbdofs;
00552    }
00553 
00554    ndofs = nvdofs + nedofs + nfdofs + nbdofs;
00555 
00556 }
00557 
00558 void FiniteElementSpace::GetElementDofs (int i, Array<int> &dofs) const
00559 {
00560    if (elem_dof)
00561    {
00562       elem_dof -> GetRow (i, dofs);
00563    }
00564    else
00565    {
00566       Array<int> V, E, Eo, F, Fo;
00567       int k, j, nv, ne, nf, nb, nfd, nd;
00568       int *ind, dim;
00569 
00570       dim = mesh->Dimension();
00571       nv = fec->DofForGeometry(Geometry::POINT);
00572       ne = (dim > 1) ? ( fec->DofForGeometry(Geometry::SEGMENT) ) : ( 0 );
00573       nb = fec->DofForGeometry(mesh->GetElementBaseGeometry(i));
00574       if (nv > 0)
00575          mesh->GetElementVertices(i, V);
00576       if (ne > 0)
00577          mesh->GetElementEdges(i, E, Eo);
00578       nfd = 0;
00579       if (dim == 3)
00580          if (fec->HasFaceDofs(mesh->GetElementBaseGeometry(i)))
00581          {
00582             mesh->GetElementFaces(i, F, Fo);
00583             for (k = 0; k < F.Size(); k++)
00584                nfd += fec->DofForGeometry(mesh->GetFaceBaseGeometry(F[k]));
00585          }
00586       nd = V.Size() * nv + E.Size() * ne + nfd + nb;
00587       dofs.SetSize(nd);
00588       if (nv > 0)
00589       {
00590          for (k = 0; k < V.Size(); k++)
00591             for (j = 0; j < nv; j++)
00592                dofs[k*nv+j] = V[k]*nv+j;
00593          nv *= V.Size();
00594       }
00595       if (ne > 0)
00596          // if (dim > 1)
00597          for (k = 0; k < E.Size(); k++)
00598          {
00599             ind = fec->DofOrderForOrientation(Geometry::SEGMENT, Eo[k]);
00600             for (j = 0; j < ne; j++)
00601                if (ind[j] < 0)
00602                   dofs[nv+k*ne+j] = -1 - ( nvdofs+E[k]*ne+(-1-ind[j]) );
00603                else
00604                   dofs[nv+k*ne+j] = nvdofs+E[k]*ne+ind[j];
00605          }
00606       ne = nv + ne * E.Size();
00607       if (nfd > 0)
00608          // if (dim == 3)
00609       {
00610          for (k = 0; k < F.Size(); k++)
00611          {
00612             ind = fec->DofOrderForOrientation(mesh->GetFaceBaseGeometry(F[k]),
00613                                               Fo[k]);
00614             nf = fec->DofForGeometry(mesh->GetFaceBaseGeometry(F[k]));
00615             for (j = 0; j < nf; j++)
00616             {
00617                if (ind[j] < 0)
00618                   dofs[ne+j] = -1 - ( nvdofs+nedofs+fdofs[F[k]]+(-1-ind[j]) );
00619                else
00620                   dofs[ne+j] = nvdofs+nedofs+fdofs[F[k]]+ind[j];
00621             }
00622             ne += nf;
00623          }
00624       }
00625       k = nvdofs + nedofs + nfdofs + bdofs[i];
00626       for (j = 0; j < nb; j++)
00627       {
00628          dofs[ne+j] = k + j;
00629       }
00630    }
00631 }
00632 
00633 const FiniteElement *FiniteElementSpace::GetFE(int i) const
00634 {
00635    const FiniteElement *FE =
00636       fec->FiniteElementForGeometry(mesh->GetElementBaseGeometry(i));
00637 
00638    if (NURBSext)
00639       NURBSext->LoadFE(i, FE);
00640 
00641    return FE;
00642 }
00643 
00644 void FiniteElementSpace::GetBdrElementDofs(int i, Array<int> &dofs) const
00645 {
00646    if (bdrElem_dof)
00647    {
00648       bdrElem_dof->GetRow(i, dofs);
00649    }
00650    else
00651    {
00652       Array<int> V, E, Eo;
00653       int k, j, nv, ne, nf, nd, iF, oF;
00654       int *ind, dim;
00655 
00656       dim = mesh->Dimension();
00657       nv = fec->DofForGeometry(Geometry::POINT);
00658       if (nv > 0)
00659          mesh->GetBdrElementVertices(i, V);
00660       ne = (dim > 1) ? ( fec->DofForGeometry(Geometry::SEGMENT) ) : ( 0 );
00661       if (ne > 0)
00662          mesh->GetBdrElementEdges(i, E, Eo);
00663       nd = V.Size() * nv + E.Size() * ne;
00664       nf = (dim == 3) ? (fec->DofForGeometry(
00665                             mesh->GetBdrElementBaseGeometry(i))) : (0);
00666       if (nf > 0)
00667       {
00668          nd += nf;
00669          mesh->GetBdrElementFace(i, &iF, &oF);
00670       }
00671       dofs.SetSize(nd);
00672       if (nv > 0)
00673       {
00674          for (k = 0; k < V.Size(); k++)
00675             for (j = 0; j < nv; j++)
00676                dofs[k*nv+j] = V[k]*nv+j;
00677          nv *= V.Size();
00678       }
00679       if (ne > 0)
00680          // if (dim > 1)
00681          for (k = 0; k < E.Size(); k++)
00682          {
00683             ind = fec->DofOrderForOrientation(Geometry::SEGMENT, Eo[k]);
00684             for (j = 0; j < ne; j++)
00685                if (ind[j] < 0)
00686                   dofs[nv+k*ne+j] = -1 - ( nvdofs+E[k]*ne+(-1-ind[j]) );
00687                else
00688                   dofs[nv+k*ne+j] = nvdofs+E[k]*ne+ind[j];
00689          }
00690       if (nf > 0)
00691          // if (dim == 3)
00692       {
00693          ne = nv + ne * E.Size();
00694          ind = (fec->DofOrderForOrientation(
00695                    mesh->GetBdrElementBaseGeometry(i), oF));
00696          for (j = 0; j < nf; j++)
00697          {
00698             if (ind[j] < 0)
00699                dofs[ne+j] = -1 - ( nvdofs+nedofs+fdofs[iF]+(-1-ind[j]) );
00700             else
00701                dofs[ne+j] = nvdofs+nedofs+fdofs[iF]+ind[j];
00702          }
00703       }
00704    }
00705 }
00706 
00707 void FiniteElementSpace::GetFaceDofs(int i, Array<int> &dofs) const
00708 {
00709    int j, k, nv, ne, nf, nd;
00710    Array<int> V, E, Eo;
00711 
00712    nv = fec->DofForGeometry(Geometry::POINT);
00713    ne = fec->DofForGeometry(Geometry::SEGMENT);
00714    if (nv > 0)
00715       mesh->GetFaceVertices(i, V);
00716    if (ne > 0)
00717       mesh->GetFaceEdges(i, E, Eo);
00718    nf = (fdofs) ? (fdofs[i+1]-fdofs[i]) : (0);
00719    nd = V.Size() * nv + E.Size() * ne + nf;
00720    dofs.SetSize (nd);
00721    if (nv > 0)
00722       for (k = 0; k < V.Size(); k++)
00723          for (j = 0; j < nv; j++)
00724             dofs[k*nv+j] = V[k]*nv+j;
00725    nv *= V.Size();
00726    // do not take edge orientation 'Eo' into account
00727    if (ne > 0)
00728       for (k = 0; k < E.Size(); k++)
00729          for (j = 0; j < ne; j++)
00730             dofs[nv+k*ne+j] = nvdofs+E[k]*ne+j;
00731    ne = nv + ne * E.Size();
00732    if (nf > 0)
00733       for (j = nvdofs+nedofs+fdofs[i], k = 0; k < nf; j++, k++)
00734          dofs[ne+k] = j;
00735 }
00736 
00737 void FiniteElementSpace::GetEdgeDofs(int i, Array<int> &dofs) const
00738 {
00739    int j, k, nv, ne;
00740    Array<int> V;
00741 
00742    nv = fec->DofForGeometry(Geometry::POINT);
00743    if (nv > 0)
00744       mesh->GetEdgeVertices(i, V);
00745    ne = fec->DofForGeometry(Geometry::SEGMENT);
00746    dofs.SetSize(2*nv+ne);
00747    if (nv > 0)
00748       for (k = 0; k < 2; k++)
00749          for (j = 0; j < nv; j++)
00750             dofs[k*nv+j] = V[k]*nv+j;
00751    nv *= 2;
00752    for (j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
00753       dofs[nv+j] = k;
00754 }
00755 
00756 void FiniteElementSpace::GetElementInteriorDofs (int i, Array<int> &dofs) const
00757 {
00758    int j, k, nb;
00759    nb = fec -> DofForGeometry (mesh -> GetElementBaseGeometry (i));
00760    dofs.SetSize (nb);
00761    k = nvdofs + nedofs + nfdofs + bdofs[i];
00762    for (j = 0; j < nb; j++)
00763    {
00764       dofs[j] = k + j;
00765    }
00766 }
00767 
00768 void FiniteElementSpace::GetEdgeInteriorDofs (int i, Array<int> &dofs) const
00769 {
00770    int j, k, ne;
00771 
00772    ne = fec -> DofForGeometry (Geometry::SEGMENT);
00773    dofs.SetSize (ne);
00774    for (j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
00775       dofs[j] = k;
00776 }
00777 
00778 const FiniteElement *FiniteElementSpace::GetBE (int i) const
00779 {
00780    const FiniteElement *BE;
00781 
00782    switch ( mesh->Dimension() )
00783    {
00784    case 1:
00785       BE = fec->FiniteElementForGeometry(Geometry::POINT);
00786    case 2:
00787       BE = fec->FiniteElementForGeometry(Geometry::SEGMENT);
00788    case 3:
00789       BE = fec->FiniteElementForGeometry(
00790          mesh->GetBdrElementBaseGeometry(i));
00791    }
00792 
00793    if (NURBSext)
00794       NURBSext->LoadBE(i, BE);
00795 
00796    return BE;
00797 }
00798 
00799 FiniteElementSpace::~FiniteElementSpace()
00800 {
00801    Destructor();
00802    // delete RefData
00803    for (int i = 0; i < RefData.Size(); i++)
00804       delete RefData[i];
00805 }
00806 
00807 void FiniteElementSpace::Destructor()
00808 {
00809    dof_elem_array.DeleteAll();
00810    dof_ldof_array.DeleteAll();
00811 
00812    if (NURBSext)
00813    {
00814       if (own_ext) delete NURBSext;
00815    }
00816    else
00817    {
00818       delete elem_dof;
00819       delete bdrElem_dof;
00820 
00821       delete [] bdofs;
00822       delete [] fdofs;
00823    }
00824 }
00825 
00826 void FiniteElementSpace::Update()
00827 {
00828    if (NURBSext)
00829    {
00830       UpdateNURBS();
00831    }
00832    else
00833    {
00834       Destructor();   // keeps RefData
00835       Constructor();
00836    }
00837 }
00838 
00839 FiniteElementSpace *FiniteElementSpace::SaveUpdate()
00840 {
00841    FiniteElementSpace *cfes = new FiniteElementSpace(*this);
00842    Constructor();
00843    return cfes;
00844 }
00845 
00846 void FiniteElementSpace::ConstructRefinementData (int k, int num_c_dofs,
00847                                                   RefinementType type)
00848 {
00849    int i,j,l;
00850    Array<int> dofs, g_dofs;
00851 
00852    RefinementData * data = new RefinementData;
00853 
00854    data -> type = type;
00855    data -> num_fine_elems = mesh -> GetNumFineElems(k);
00856 
00857    // we assume that each fine element has <= dofs than the
00858    // initial coarse element
00859    data -> fl_to_fc = new Table(data -> num_fine_elems, num_c_dofs);
00860    for (i = 0; i < data -> num_fine_elems; i++) {
00861       GetElementDofs(mesh -> GetFineElem(k,i), g_dofs); // TWO_LEVEL_FINE
00862       for (j = 0; j < g_dofs.Size(); j++)
00863          data -> fl_to_fc -> Push (i,dofs.Union(g_dofs[j]));
00864    }
00865    data -> fl_to_fc -> Finalize();
00866    data -> num_fine_dofs = dofs.Size();
00867 
00868    // construction of I
00869 
00870    // k is a coarse element index but mesh is in TWO_LEVEL_FINE state
00871    int geomtype = mesh->GetElementBaseGeometry(k);
00872    const FiniteElement *fe = fec -> FiniteElementForGeometry(geomtype);
00873    // const IntegrationRule &ir = fe -> GetNodes();
00874    int nedofs = fe -> GetDof();  // number of dofs for each element
00875 
00876    ElementTransformation *trans;
00877    // IntegrationPoint ip;
00878    // DenseMatrix tr;
00879    Array<int> row;
00880 
00881    // Vector shape(nedofs);
00882    DenseMatrix I (nedofs);
00883    data -> I = new DenseMatrix(data -> num_fine_dofs, nedofs);
00884 
00885    for (i=0; i < data -> num_fine_elems; i++) {
00886 
00887       trans = mesh -> GetFineElemTrans(k,i);
00888       // trans -> Transform(ir,tr);
00889       fe -> GetLocalInterpolation (*trans, I);
00890       data -> fl_to_fc -> GetRow(i,row);
00891 
00892       for (j=0; j < nedofs; j++)
00893       {
00894          /*
00895            ip.x = tr(0,j); ip.y = tr(1,j);       if (tr.Height()==3) ip.z = tr(2,j);
00896            fe -> SetCalcElemTrans(trans); shape(0) = j;
00897            fe -> CalcShape(ip, shape);
00898          */
00899          for (l=0; l < nedofs; l++)
00900             (*(data->I))(row[j],l) = I(j,l);
00901       }
00902    }
00903 
00904    RefData.Append(data);
00905 }
00906 
00907 void FiniteElementSpace::Save(ostream &out) const
00908 {
00909    out << "FiniteElementSpace\n"
00910        << "FiniteElementCollection: " << fec->Name() << '\n'
00911        << "VDim: " << vdim << '\n'
00912        << "Ordering: " << ordering << '\n';
00913 }
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines