00001 #ifndef _MY_DealBaseDoFs_
00002 #define _MY_DealBaseDoFs_
00003 #include "dealbase.h"
00004 #include <dofs/dof_handler.h>
00005 #include "globals.h"
00006
00007
00008
00009 class DealBaseDoFs
00010 {
00011 private:
00012
00013 protected:
00014
00015 public:
00016 DealBaseDoFs(){}
00017 ~DealBaseDoFs(){}
00018
00019 static DoFCell getDoFCellAtPoint(const DoFHandler<3>& dof,const Point<3> &p);
00020
00021 static DoFCell getAdjDoFCell(DoFCell cell,CellDirection3D dir);
00022
00023 static double getDoFVerticeValue(DoFCell cell,VertexDirection3D dir,unsigned deg,const VecDouble &vSol);
00024 void projectDofSolutionAtVertice(const DoFHandler<3>& dofs,int deg,const VecDouble &vDofSol,VecDouble &vSol);
00025 void printCellDofsInformation(const DoFHandler<2>& dofs,std::ostream &out);
00026 void printCellDofsInformation(const DoFHandler<3>& dofs);
00027 static void getVerticeValuesFromDofs(DoFHandler<3> &dofs,const VecDouble &vSol,VecDouble &vVerticesValues,unsigned deg);
00028
00029 void getCentralValuesFromDofs(DoFHandler<3> &dof,const VecDouble &sol, VecDouble &cValues,unsigned component);
00030
00031 void printFE(FiniteElement<3> &fe);
00032 void getValuesAtCells(DoFHandler<3> &dofs,const VecDouble &sol, VecDouble cValues);
00033
00034
00035 template<class Vetor>
00036 void getValuesAtFaces(DoFHandler<3> &dof,const Vetor &sol, Matrix &Mvalues,const std::vector<unsigned> &components, const bool homMesh);
00037
00038
00039
00040
00041
00042 };
00043
00044
00045
00046
00064 template<class Vetor>
00065 void DealBaseDoFs::getValuesAtFaces(DoFHandler<3> &dof,const Vetor &sol, Matrix &Mvalues,const std::vector<unsigned> &components, const bool homMesh)
00066 {
00067 assert(Mvalues.n() == components.size());
00068 assert(Mvalues.m() == dof.get_tria().n_active_faces());
00069 std::vector<Point<3> > points;
00070 points.push_back(Point<3>(0,0.5,0.5));
00071 points.push_back(Point<3>(1,0.5,0.5));
00072 points.push_back(Point<3>(0.5,0,0.5));
00073 points.push_back(Point<3>(0.5,1,0.5));
00074 points.push_back(Point<3>(0.5,0.5,0));
00075 points.push_back(Point<3>(0.5,0.5,1));
00076 std::vector<double> weights(6,1);
00077 unsigned n_quadrature_points = points.size();
00078
00079 Quadrature<3> quad(points,weights);
00080
00081 FEValues<DIM> fe_values(dof.get_fe(),quad,update_values);
00082
00083 const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
00084 const unsigned int n_components = components.size();
00085 std::vector<unsigned int> local_dof_indices (dofs_per_cell);
00086
00087
00088 DoFHandler<3>::active_cell_iterator cell = dof.begin_active(),
00089 endc = dof.end();
00090
00091
00092 if (homMesh)
00093 fe_values.reinit (cell);
00094 for (; cell!=endc; ++cell)
00095 {
00096 if (!homMesh)
00097 fe_values.reinit (cell);
00098
00099 cell->get_dof_indices(local_dof_indices);
00100
00101
00102 for (unsigned q_point=0;q_point < n_quadrature_points;q_point++)
00103 {
00104
00105
00106
00107 double *values = &(Mvalues(cell->face_index(q_point),0));
00108
00109
00110 for (unsigned j=0;j < n_components;j++)
00111 {
00112 double acum = 0.0;
00113
00114 for (unsigned int i=0; i<dofs_per_cell; ++i)
00115 {
00116 acum += sol(local_dof_indices[i])*fe_values.shape_value_component(i,q_point,components[j]);
00117 }
00118 values[j]=acum;
00119
00120 }
00121 }
00122 }
00123 }
00124
00125
00126 #endif