MFEM v2.0
|
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 }