00001 #include "dealbase.h"
00002 #include "globals.h"
00003 #include "numericmethods.h"
00004
00005 const int DealBase::INVALID_INDEX = -1;
00006
00007
00014 void DealBase::setTriangulation(Triangulation<DIM>& trig)
00015 {
00016 m_pTriangulation = &trig;
00017 END_CELL = getTriangulation().end();
00018 n_cells=getTriangulation().n_active_cells();
00019 }
00020
00021
00022
00029 double DealBase::getCentralValue(const Cell &cell,const VecDouble &values) const
00030 {
00031 assert(values.size() == numCells());
00032 return values(cell->index());
00033 }
00034
00035
00045 double DealBase::getVerticeValue(const Cell &cell,VertexDirection3D vertex,const VecDouble &values)
00046 {
00047 assert(values.size() == numVertices());
00048 return values(cell->vertex_index(vertex));
00049 }
00050
00051
00059 unsigned DealBase::getCentralPointIndex(Cell &cell)
00060 {
00061 return cell->index();
00062 }
00063
00064
00065
00066 void DealBase::setCentralValue(const Cell &cell,double dd,VecDouble &values)
00067 {
00068 values(cell->index()) = dd;
00069 }
00070
00071 void DealBase::setCentralValue(DoFCell &cell,double dd,VecDouble &values)
00072 {
00073 values(cell->index()) = dd;
00074 }
00075
00076
00077 void DealBase::setVerticesValuesFromFunction(Function<DIM> &f,VecDouble &vValues,unsigned deg)
00078 {
00079 assert(vValues.size() == numVertices());
00080 Cell cell = getTriangulation().begin_active(),
00081 endc = getTriangulation().end();
00082 for (;cell!=endc;++cell)
00083 {
00084 for (unsigned i=0;i<GeometryInfo<DIM>::vertices_per_cell;i++)
00085 {
00086 setVerticeValue(cell,(VertexDirection3D) i,f.value(cell->vertex(i),deg),vValues);
00087 }
00088 }
00089 }
00090
00091
00094 void DealBase::getVerticePoint(unsigned vIndex,Point<DIM>& X) const
00095 {
00096 assert(vIndex < getTriangulation().n_vertices());
00097 X = getTriangulation().get_vertices()[vIndex];
00098 }
00099
00100 void DealBase::printVertices()
00101 {
00102 unsigned i;
00103 for (i=0;i<getTriangulation().n_vertices();i++)
00104 printf("%d) (%g,%g)\n",i,getTriangulation().get_vertices()[i](0),getTriangulation().get_vertices()[i](1));
00105 }
00106
00107 void DealBase::printCells()
00108 {
00109 Cell cell = getTriangulation().begin_active();
00110 Cell endc = getTriangulation().end();
00111
00112 for (;cell!=endc;cell++)
00113 {
00114 printf("Cell %d) CC=(%g,%g), VInd: %d %d %d %d\n",cell->index(),cell->barycenter()(0),cell->barycenter()(1),cell->vertex_index(0),cell->vertex_index(1),cell->vertex_index(2),cell->vertex_index(3));
00115
00116 }
00117 }
00118
00119
00120
00121 void DealBase::printMatrixMaple(SparseMatrix<double> &M)
00122 {
00123 unsigned i,j;
00124
00125 cout << "Matrix("<<M.m()<<","<<M.n() << ",[";
00126 for (i=0;i<M.m();i++)
00127 {
00128 cout << "[";
00129 for (j=0;j<M.n()-1;j++)
00130 {
00131 cout << M.el(i,j) << ",";
00132 }
00133 if (i == M.m()-1)
00134 cout << M.el(i,j) << "]])";
00135 else
00136 cout << M.el(i,j) << "],";
00137 }
00138 cout.flush();
00139
00140 }
00141
00142
00143 void DealBase::printMatrixMaple(FullMatrix<double> &M)
00144 {
00145 unsigned i,j;
00146
00147 cout << "Matrix("<<M.m()<<","<<M.n() << ",[";
00148 for (i=0;i<M.m();i++)
00149 {
00150 cout << "[";
00151 for (j=0;j<M.n()-1;j++)
00152 {
00153 cout << M(i,j) << ",";
00154 }
00155 if (i == M.m()-1)
00156 cout << M(i,j) << "]])";
00157 else
00158 cout << M(i,j) << "],";
00159 }
00160 cout.flush();
00161
00162 }
00163
00164
00165
00172 void DealBase::printPoint(const Point<DIM> &P,std::ostream &out)
00173 {
00174 out << "[ ";
00175 unsigned i;
00176 for (i=0;i<DIM;i++)
00177 out << P(i) << " ";
00178 out << "]";
00179 }
00180
00181
00182
00190 void DealBase::setMaterialIdFromFunctionDomain(GeneralFunctionInterface &f,unsigned material_id)
00191 {
00192 Cell cell = getTriangulation().begin_active();
00193 Cell endc = getTriangulation().end();
00194
00195 for (;cell!=endc;cell++)
00196 {
00197 if (f.isInDomain(cell->barycenter()))
00198 {
00199 cell->set_material_id(material_id);
00200 cout << "done \n";
00201 }
00202
00203 }
00204 }
00205
00206
00207
00208
00216 void DealBase::setMaterialIdOfPrescribedBoundaryCells(GeneralFunctionInterface &f,unsigned material_id)
00217 {
00218 Triangulation<DIM> &tria = getTriangulation();
00219 Cell cell = tria.begin_active();
00220 Cell endc = tria.end();
00221
00222
00223 for (;cell!=endc;cell++)
00224 {
00225 if (cell->at_boundary())
00226 {
00227 for (unsigned i=0;i<GeometryInfo<DIM>::vertices_per_cell;i++)
00228 {
00229 if (f.isInDomain(cell->vertex(i)))
00230 {
00231 cell->set_material_id(material_id);
00232 break;
00233 }
00234 }
00235 }
00236 }
00237 }
00238
00239
00240
00241
00248 void DealBase::setVerticeValue(Cell &cell,VertexDirection3D vertex,double dd,VecDouble &values)
00249 {
00250 assert(values.size() == numVertices());
00251 values(cell->vertex_index(vertex)) = dd;
00252
00253 }
00254
00255
00256
00257
00266 void DealBase::setPrescribedBoundaryId(Function3D &f,unsigned boundary_id)
00267 {
00268 Cell cell = getTriangulation().begin_active();
00269 Cell endc = getTriangulation().end();
00270
00271 for (;cell!=endc;cell++)
00272 {
00273 for (unsigned face_no=0;face_no < GeometryInfo<DIM>::faces_per_cell;face_no++)
00274 {
00275 Face face = cell->face(face_no);
00276 if (face->at_boundary())
00277 {
00278 bool markFace = true;
00279 for (unsigned ivertex=0;ivertex< GeometryInfo<DIM>::vertices_per_face;ivertex++)
00280 {
00281 if (!f.isInDomain(face->vertex(ivertex)))
00282 {
00283 markFace = false;
00284 break;
00285 }
00286 }
00287 if (markFace)
00288 {
00289 face->set_boundary_indicator(boundary_id);
00290 }
00291 }
00292 }
00293 }
00294 }
00295
00296
00297 void DealBase::setCellFlag(const Cell &cell,const bool flag,VecBool &cVFlags) const
00298 {
00299 assert(cVFlags.size() == numCells());
00300 cVFlags[cell->index()] = flag;
00301 }
00302
00303 bool DealBase::getCellFlag(const Cell &cell,const VecBool &cVFlags) const
00304 {
00305 assert(cVFlags.size() == numCells());
00306 return cVFlags[cell->index()];
00307
00308 }
00309
00310
00311 void DealBase::setVerticeFlag(const Cell &cell,VertexDirection3D dir,const bool flag,VecBool &vVFlags) const
00312 {
00313 assert(vVFlags.size() == numVertices());
00314 vVFlags[cell->vertex_index(dir)] = flag;
00315 }
00316
00317 bool DealBase::getVerticeFlag(const Cell &cell,VertexDirection3D dir,const VecBool &vVFlags) const
00318 {
00319 assert(vVFlags.size() == numVertices());
00320 return vVFlags[cell->vertex_index(dir)];
00321 }
00322
00323 void DealBase::printVector(const VecDouble &vSol)
00324 {
00325 for (unsigned i=0;i<vSol.size();i++)
00326 printf("%d) %g\n",i,vSol(i));
00327 }
00328
00329
00330
00339 double DealBase::getNeighborValue(const Cell &cell,CellDirection3D dir,const VecDouble &cValues)
00340 {
00341 assert(cell->neighbor_index(dir)!=INVALID_INDEX);
00342 assert(cValues.size() == numCells());
00343 return cValues(cell->neighbor_index(dir));
00344 }
00345
00346
00354 Cell DealBase::getCellAtPoint(const Point<DIM> &p)
00355 {
00356 Cell cell = beginCell();
00357 Cell endc = endCell();
00358 for (;cell!= endc;cell++)
00359 {
00360 if (NumericMethods::isInCube(p,cell->vertex(VERTEX_000),cell->vertex(VERTEX_111)))
00361 {
00362 return cell;
00363 }
00364 }
00365 return Cell();
00366
00367 }
00368
00369
00370 double DealBase::getFaceValue(const Face &face,const VecDouble &values)
00371 {
00372 assert(values.size() == numFaces());
00373 return values(face->index());
00374 }
00375
00376
00377 void DealBase::setFaceValue(const Face &face,double dd,VecDouble &values)
00378 {
00379 assert(values.size() == numFaces());
00380 values(face->index())=dd;
00381 }
00382
00383 double DealBase::getFaceValue(const Cell &cell,FaceDirection3D dir,const VecDouble &values)
00384 {
00385 assert(values.size() == numFaces());
00386 return values(cell->face(dir)->index());
00387 }
00388
00389 unsigned DealBase::getFaceValueIndex(const Cell &cell,FaceDirection3D dir)
00390 {
00391 return cell->face(dir)->index();
00392 }
00393
00394 void DealBase::setFaceValue(const Cell &cell,FaceDirection3D dir,double dd,VecDouble &values)
00395 {
00396 assert(values.size() == numFaces());
00397 values(cell->face(dir)->index())=dd;
00398 }
00399
00400
00401
00404 void DealBase::printFaces()
00405 {
00406 Face face = getTriangulation().begin_active_face();
00407 Face endFace = getTriangulation().end_face();
00408
00409 for(;face!=endFace;face++)
00410 {
00411 printf("Face %d) Vertices Indices: %d,%d\n",face->index(),face->vertex_index(0),face->vertex_index(1));
00412 }
00413 }
00414
00415 void DealBase::setFaceFlag(const Face &face,bool b,VecBool &values)
00416 {
00417 assert(values.size() == numFaces());
00418 values[face->index()]=b;
00419
00420 }
00421 bool DealBase::getFaceFlag(const Face &face,const VecBool& values)
00422 {
00423 assert(values.size() == numFaces());
00424 return values[face->index()];
00425 }
00426
00427 void DealBase::flagPrescribedFaces(Function3D &f,VecBool &vFlags)
00428 {
00429 assert(vFlags.size() == numFaces());
00430 for (Face face=beginFace();face!=endFace();face++)
00431 {
00432 Point3D p = face->center();
00433 setFaceFlag(face,face->at_boundary() && (f.isInDomain(p)),vFlags);
00434 }
00435 }
00436
00437 void DealBase::printPoint(const Point<DIM> &p)
00438 {
00439 if (DIM == 2)
00440 printf("<%5g,%5g>",p[0],p[1]);
00441 else if (DIM == 1)
00442 printf("<%5g>",p[0]);
00443 else if (DIM == 3)
00444 printf("<%5g,%5g,%5g>",p[0],p[1],p[2]);
00445 else
00446 abort();
00447 }
00448
00449 void DealBase::printCentralValues(const VecDouble& cValues)
00450 {
00451 assert(cValues.size()==numCells());
00452 Cell cell = beginCell();
00453 Cell endc = endCell();
00454
00455 for (;cell!=endc;cell++)
00456 {
00457 printf("Cell: <%3g,%3g> %5g\n",cell->center()(0),cell->center()(1),getCentralValue(cell,cValues));
00458 }
00459 }
00460 void DealBase::printCentralValues(const VecDouble& cValues1,const VecDouble& cValues2)
00461 {
00462 assert(cValues1.size()==numCells());
00463 assert(cValues2.size()==numCells());
00464 Cell cell = beginCell();
00465 Cell endc = endCell();
00466
00467 for (;cell!=endc;cell++)
00468 {
00469 printf("Cell: <%3g,%3g> %5g %5g\n",cell->center()(0),cell->center()(1),getCentralValue(cell,cValues1),getCentralValue(cell,cValues2));
00470 }
00471 }
00472
00473 void DealBase::printCentralValues(const VecDouble& cValues1,const VecDouble& cValues2,const VecDouble& cValues3,const VecDouble& cValues4)
00474 {
00475 assert(cValues1.size()==numCells());
00476 assert(cValues2.size()==numCells());
00477 assert(cValues3.size()==numCells());
00478 assert(cValues4.size()==numCells());
00479 Cell cell = beginCell();
00480 Cell endc = endCell();
00481
00482 for (;cell!=endc;cell++)
00483 {
00484 printf("Cell: <%3g,%3g> %5g %5g %5g %5g\n",cell->center()(0),cell->center()(1),getCentralValue(cell,cValues1),getCentralValue(cell,cValues2),getCentralValue(cell,cValues3),getCentralValue(cell,cValues4));
00485 }
00486 }
00487
00488
00489
00490
00491 void DealBase::printFacesValuesAtBoundary(const VecDouble& fValues)
00492 {
00493 assert(fValues.size()==numFaces());
00494 Cell cell = beginCell();
00495 Cell endc = endCell();
00496
00497 for (;cell!=endc;cell++)
00498 {
00499 for (unsigned face_no=0;face_no<GeometryInfo<DIM>::faces_per_cell;face_no++)
00500 {
00501 Face face = cell->face(face_no);
00502 if (face->at_boundary())
00503 printf("Face: <%3g,%3g> %5g\n",face->center()(0),face->center()(1),getFaceValue(face,fValues));
00504 }
00505 }
00506 }
00507
00508
00509
00510 void DealBase::printFacesValuesAtBoundary(const VecDouble& fValues,const VecDouble& fValues2)
00511 {
00512 assert(fValues.size()==numFaces());
00513 assert(fValues2.size()==numFaces());
00514 Cell cell = beginCell();
00515 Cell endc = endCell();
00516
00517 for (;cell!=endc;cell++)
00518 {
00519 for (unsigned face_no=0;face_no<GeometryInfo<DIM>::faces_per_cell;face_no++)
00520 {
00521 Face face = cell->face(face_no);
00522 if (face->at_boundary())
00523 printf("Face: <%3g,%3g> %5g %5g\n",face->center()(0),face->center()(1),getFaceValue(face,fValues),getFaceValue(face,fValues2));
00524 }
00525 }
00526 }
00527
00528
00529 void DealBase::printFacesValuesAtBoundary(const VecDouble& fValues1,const VecDouble& fValues2,const VecDouble& fValues3,const VecDouble& fValues4)
00530 {
00531 assert(fValues1.size()==numFaces());
00532 assert(fValues2.size()==numFaces());
00533 assert(fValues3.size()==numFaces());
00534 assert(fValues4.size()==numFaces());
00535 Cell cell = beginCell();
00536 Cell endc = endCell();
00537
00538 for (;cell!=endc;cell++)
00539 {
00540 for (unsigned face_no=0;face_no<GeometryInfo<DIM>::faces_per_cell;face_no++)
00541 {
00542 Face face = cell->face(face_no);
00543 if (face->at_boundary())
00544 printf("Face: <%3g,%3g> %5g %5g %5g %5g\n",face->center()(0),face->center()(1),getFaceValue(face,fValues1),getFaceValue(face,fValues2),getFaceValue(face,fValues3),getFaceValue(face,fValues4));
00545 }
00546 }
00547
00548 }
00549
00550
00551
00559 Cell DealBase::getAdjCell(Cell &cell,CellDirection3D dir)
00560 {
00561 if (cell->face(dir)->at_boundary())
00562 return getEndCell();
00563 else
00564 return cell->neighbor(dir);
00565 }
00566
00567
00568
00578 Cell DealBase::getAdjCell(Cell &cell,CellDirection3D dir1,CellDirection3D dir2)
00579 {
00580 static Cell auxCell;
00581 if (cell->face(dir1)->at_boundary())
00582 return getEndCell();
00583 else
00584 {
00585 auxCell = cell->neighbor(dir1);
00586 if (auxCell->face(dir2)->at_boundary())
00587 return getEndCell();
00588 else
00589 return auxCell->neighbor(dir2);
00590 }
00591 }
00592
00593
00594
00595
00596
00597 void DealBase::printCell(Cell cell,std::ostream &out)
00598 {
00599 char str[500];
00600 sprintf(str,"==== Cell %d ====\n",cell->index());
00601 out << str;
00602 out << "Barycenter: " << cell->barycenter() << endl;
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612 out << "V000 = " << cell->vertex_index(VERTEX_000) << ") " << cell->vertex(VERTEX_000) << endl;
00613 out << "V100 = " << cell->vertex_index(VERTEX_100) << ") " << cell->vertex(VERTEX_100) << endl;
00614 out << "V110 = " << cell->vertex_index(VERTEX_110) << ") " << cell->vertex(VERTEX_110) << endl;
00615 out << "V010 = " << cell->vertex_index(VERTEX_010) << ") " << cell->vertex(VERTEX_010) << endl;
00616 out << "V001 = " << cell->vertex_index(VERTEX_001) << ") " << cell->vertex(VERTEX_001) << endl;
00617 out << "V101 = " << cell->vertex_index(VERTEX_101) << ") " << cell->vertex(VERTEX_101) << endl;
00618 out << "V111 = " << cell->vertex_index(VERTEX_111) << ") " << cell->vertex(VERTEX_111) << endl;
00619 out << "V011 = " << cell->vertex_index(VERTEX_011) << ") " << cell->vertex(VERTEX_011) << endl;
00620
00621 out << "RIGHT_FACE " << cell->face_index (RIGHT_FACE) << ")" << endl;
00622 out << "LEFT_FACE " << cell->face_index (LEFT_FACE) << ")" << endl;
00623 out << "FRONT_FACE " << cell->face_index (FRONT_FACE) << ")" << endl;
00624 out << "BACK_FACE " << cell->face_index (BACK_FACE) << ")" << endl;
00625 out << "BOTTOM_FACE " << cell->face_index(BOTTOM_FACE) << ")"<< endl;
00626 out << "UP_FACE " << cell->face_index (UP_FACE) << ")" << endl;
00627
00628
00629 }
00630
00631 void DealBase::printTriangulation(Triangulation<3> &tria,std::ostream &out)
00632 {
00633 Cell cell = tria.begin_active();
00634 Cell endc = tria.end();
00635 char str[500];
00636 for (;cell != endc;cell++)
00637 {
00638 sprintf(str,"==== Cell %d ====\n",cell->index());
00639 out << str;
00640 out << "Barycenter: " << cell->barycenter() << endl;
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650 out << "V000 = " << cell->vertex_index(VERTEX_000) << ") " << cell->vertex(VERTEX_000) << endl;
00651 out << "V100 = " << cell->vertex_index(VERTEX_100) << ") " << cell->vertex(VERTEX_100) << endl;
00652 out << "V110 = " << cell->vertex_index(VERTEX_110) << ") " << cell->vertex(VERTEX_110) << endl;
00653 out << "V010 = " << cell->vertex_index(VERTEX_010) << ") " << cell->vertex(VERTEX_010) << endl;
00654 out << "V001 = " << cell->vertex_index(VERTEX_001) << ") " << cell->vertex(VERTEX_001) << endl;
00655 out << "V101 = " << cell->vertex_index(VERTEX_101) << ") " << cell->vertex(VERTEX_101) << endl;
00656 out << "V111 = " << cell->vertex_index(VERTEX_111) << ") " << cell->vertex(VERTEX_111) << endl;
00657 out << "V011 = " << cell->vertex_index(VERTEX_011) << ") " << cell->vertex(VERTEX_011) << endl;
00658
00659 out << "RIGHT_FACE " << cell->face_index (RIGHT_FACE) << ")" << endl;
00660 out << "LEFT_FACE " << cell->face_index (LEFT_FACE) << ")" << endl;
00661 out << "FRONT_FACE " << cell->face_index (FRONT_FACE) << ")" << endl;
00662 out << "BACK_FACE " << cell->face_index (BACK_FACE) << ")" << endl;
00663 out << "BOTTOM_FACE " << cell->face_index(BOTTOM_FACE) << ")"<< endl;
00664 out << "UP_FACE " << cell->face_index (UP_FACE) << ")" << endl;
00665
00666 }
00667 }
00668
00675 int DealBase::getAdjCellIndex(unsigned iCell,CellDirection3D dir)
00676 {
00677 assert(iCell < numCells());
00678 return CellAccessor<3>(&getTriangulation(),0,iCell).neighbor_index(dir);
00679 }
00680
00681
00689 double DealBase::getAdjCellValue(unsigned iCell,CellDirection3D dir,const VecDouble &v)
00690 {
00691 assert(getAdjCellIndex(iCell,dir) !=INVALID_INDEX);
00692 assert(v.size() == this->numCells());
00693 return v(getAdjCellIndex(iCell,dir));
00694 }
00695
00696
00697
00698
00699
00704 bool DealBase::isCellAtBoundary(unsigned iCell)
00705 {
00706 return CellAccessor<3>(&getTriangulation(),0,iCell).at_boundary();
00707 }
00708
00709
00710
00716 bool DealBase::hasAdjCell(unsigned iCell, CellDirection3D dir)
00717 {
00718 return getAdjCellIndex(iCell,dir) == INVALID_INDEX;
00719 }
00720
00726 int DealBase::getFaceIndex(unsigned iCell,FaceDirection3D dir)
00727 {
00728 return CellAccessor<3>(&getTriangulation(),0,iCell).quad_index(dir);
00729
00730 }
00731
00732
00733
00740 double DealBase::getFaceValue(unsigned iCell,FaceDirection3D dir, const VecDouble &v)
00741 {
00742 assert(getFaceIndex(iCell,dir) != INVALID_INDEX);
00743 assert(v.size() == numFaces());
00744 return v(getFaceIndex(iCell,dir));
00745 }
00746
00753 void DealBase::setCentralValuesFromFunction(Function3D &f,VecDouble &cValues,unsigned component)
00754 {
00755 assert(cValues.size() == numCells());
00756
00757 Cell cell = getTriangulation().begin_active(),
00758 endc = getTriangulation().end();
00759 Point<DIM> BC;
00760 for (;cell!=endc;++cell)
00761 {
00762 BC = cell->barycenter();
00763 setCentralValue(cell,f(BC,component),cValues);
00764 }
00765 }
00773 void DealBase::setCentralValuesInFunctionDomain(Function3D &f,MapIntDouble &values,unsigned component)
00774 {
00775
00776 Cell cell = getTriangulation().begin_active(),
00777 endc = getTriangulation().end();
00778 Point<DIM> BC;
00779 for (;cell!=endc;++cell)
00780 {
00781 BC = cell->barycenter();
00782 if (f.isInDomain(BC,component))
00783 values[cell->index()] = f(BC,component);
00784 }
00785 }
00786
00787
00788
00797 void DealBase::setFacesValuesFromFunction(Function3D &f,VecDouble &fValues,unsigned deg,bool bOnlyPrescribedBoundary)
00798 {
00799 assert(fValues.size() == numFaces());
00800 Face face = getTriangulation().begin_active_face();
00801 Face endFace = getTriangulation().end_face();
00802
00803 for (;face!=endFace;face++)
00804 {
00805 if (bOnlyPrescribedBoundary && !face->at_boundary())
00806 continue;
00807
00808 Point<3> P = face->vertex(BL_VERTEX);
00809 P+=face->vertex(UR_VERTEX);
00810 P[0]=P[0]/2.0;
00811 P[1]=P[1]/2.0;
00812 P[2]=P[2]/2.0;
00813 if (!f.isInDomain(P,deg))
00814 continue;
00815 else
00816 setFaceValue(face,f(P,deg),fValues);
00817 }
00818 }
00819
00820
00821
00828 void DealBase::tagFacesInDomain(Function3D &f,VecTag &tags,char value)
00829 {
00830 assert(tags.size() == numFaces());
00831 Face face = getTriangulation().begin_active_face();
00832 Face endFace = getTriangulation().end_face();
00833
00834 for (;face!=endFace;face++)
00835 {
00836 Point<3> P = face->vertex(BL_VERTEX);
00837 P+=face->vertex(UR_VERTEX);
00838 P[0]=P[0]/2.0;
00839 P[1]=P[1]/2.0;
00840 P[2]=P[2]/2.0;
00841 if (f.isInDomain(P,0))
00842 tags[face->index()] = value;
00843 }
00844 }
00845
00846
00847 void DealBase::tagFacesInDomain(Function3D &f,VecBool &tags)
00848 {
00849 assert(tags.size() == numFaces());
00850 Face face = getTriangulation().begin_active_face();
00851 Face endFace = getTriangulation().end_face();
00852
00853 for (;face!=endFace;face++)
00854 {
00855 Point<3> P = face->vertex(BL_VERTEX);
00856 P+=face->vertex(UR_VERTEX);
00857 P[0]=P[0]/2.0;
00858 P[1]=P[1]/2.0;
00859 P[2]=P[2]/2.0;
00860 if (f.isInDomain(P,0))
00861 tags[face->index()] = true;
00862 }
00863 }
00864
00865
00866
00867
00875 void DealBase::setFacesValuesFromFunction(Function3D &f,Matrix &MValues,bool bOnlyAtBoundary)
00876 {
00877 assert(MValues.m() == numFaces());
00878 assert(MValues.n() == f.n_components());
00879
00880 Face face = getTriangulation().begin_active_face();
00881 Face endFace = getTriangulation().end_face();
00882
00883 for (;face!=endFace;face++)
00884 {
00885 if (bOnlyAtBoundary && !face->at_boundary())
00886 continue;
00887
00888 Point<3> P = face->vertex(BL_VERTEX);
00889 P+=face->vertex(UR_VERTEX);
00890 P[0]=P[0]/2.0;
00891 P[1]=P[1]/2.0;
00892 P[2]=P[2]/2.0;
00893 unsigned faceIndex = face->index();
00894 for (unsigned deg=0;deg<f.n_components();deg++)
00895 {
00896 MValues(faceIndex,deg) = f(P,deg);
00897 }
00898 }
00899 }
00900
00901
00902 double DealBase::faceArea(Face face)
00903 {
00904 return NumericMethods::squareMeasure(face->vertex(BL_VERTEX),
00905 face->vertex(BR_VERTEX),
00906 face->vertex(UR_VERTEX),
00907 face->vertex(UL_VERTEX));
00908 }
00909
00910
00911
00912 void DealBase::getFaceBarycenter(Face &face,Point3D &p)
00913 {
00914 p[0]=0;
00915 p[1]=0;
00916 p[2]=0;
00917
00918 p+=face->vertex(0);
00919 p+=face->vertex(1);
00920 p+=face->vertex(2);
00921 p+=face->vertex(3);
00922 p/=4.0;
00923 }
00924
00925 Point3D DealBase::getFaceBarycenter(Face &face)
00926 {
00927 Point3D p;
00928 getFaceBarycenter(face,p);
00929 return p;
00930 }
00931
00932
00938 void DealBase::getCellsInFunctionDomain(Function3D &f,VecIndex &vec)
00939 {
00940 vec.clear();
00941 Cell cell = getTriangulation().begin_active(),
00942 endc = getTriangulation().end();
00943 for (;cell!=endc;++cell)
00944 {
00945 Point3D barycenter = cell->barycenter();
00946 if (f.isInDomain(barycenter))
00947 {
00948 vec.push_back(cell->index());
00949 }
00950
00951 }
00952 }
00953
00954
00955
00956
00957
00958
00959
00960 Cell DealBase::getCell(int index)
00961 {
00962 return Cell(&getTriangulation(),beginCell()->level(),index);
00963 }
00964
00965
00966 void DealBase::printCells(VecIndex &vec,std::ostream &out)
00967 {
00968 for (unsigned i = 0;i<vec.size();i++)
00969 {
00970 out << i << ")";
00971 printCell(getCell(vec[i]),out);
00972 }
00973 }
