00001 #include "varformpoisson.h"
00002 #include "base/quadrature_lib.h"
00003 #include "unitcube.h"
00004 VarFormPoisson::VarFormPoisson( FEOrthoMesh &fe)
00005 :VarFormOrthoMesh(fe),_feValues(fe)
00006 {
00007 assert(_fe.n_fields()==1);
00008 _MSparsity.reinit(_fe.n_local_dofs(),_fe.n_local_dofs());
00009 double c=1.0;
00010 _MSparsity.fill(c);
00011 dealii::QGauss<3> quad(2);
00012 dealii::QGauss<2> quad2d(2);
00013
00014 _feValues.setPoints(quad);
00015 pvK=NULL;
00016 _vFaceValues.resize(UnitCube::n_faces());
00017
00018
00019 for (unsigned face_dir=0;face_dir<_vFaceValues.size();face_dir++)
00020 {
00021 _vFaceValues[face_dir]=new MappedFEFaceValues(fe,face_dir);
00022 _vFaceValues[face_dir]->setQuadrature(quad2d);
00023 }
00024 }
00025
00026 void VarFormPoisson::assembly_local_system(OrthoMesh::Cell_It &cell,Matrix &ML)
00027 {
00028 assert(pvK);
00029 assert(pvK->size() == getMesh()->numCells());
00030 unsigned n_local_dofs = _fe.n_local_dofs();
00031 unsigned n_pts = _feValues.n_pts();
00032 ML=0.0;
00033 double K = (*pvK)(cell->index());
00034 for (unsigned i=0;i<n_local_dofs;i++)
00035 {
00036 for (unsigned j=0;j<n_local_dofs;j++)
00037 {
00038 for (unsigned ptId=0;ptId<n_pts;ptId++)
00039 {
00040 ML(i,j)+=_feValues.function_grad(ptId,i)*_feValues.function_grad(ptId,j)*_feValues.JxW(ptId);
00041 }
00042 }
00043 }
00044 ML*=K;
00045 }
00046
00047 void VarFormPoisson::assembly_local_rhs(OrthoMesh::Cell_It &cell,VecDouble &BL)
00048 {
00049 BL=0.0;
00050 }
00051
00052
00053
00054 void VarFormPoisson::mesh_change(OrthoMesh &mesh)
00055 {
00056 _map.reinit(mesh);
00057 _feValues.compute_values();
00058 _feValues.compute_grads(_map);
00059 _feValues.compute_JxW(_map);
00060
00061 for (unsigned i=0;i<UnitCube::n_faces();i++)
00062 {
00063 _vFaceValues[i]->compute_values();
00064
00065 _vFaceValues[i]->compute_JxW(_map);
00066 }
00067
00068 }
00069
00070
00071
00072 void VarFormPoisson::neumann_condition(VecDouble &Bl,OrthoMesh::Cell_It cell,OrthoMesh::Face_It face,unsigned face_dir,Function3D &fN)
00073 {
00074
00075 MappedFEFaceValues &fV = *( _vFaceValues[face_dir]);
00076 _map.reinit(cell);
00077 fV.compute_points(_map);
00078
00079 unsigned n_local_dofs = _fe.face_to_cell_local_map(face_dir).size();
00080 unsigned n_pts = fV.n_pts();
00081 Bl=0.0;
00082 for (unsigned i=0;i<n_local_dofs;i++)
00083 for (unsigned ptId=0;ptId<n_pts;ptId++)
00084 {
00085 Bl(i)+=fV.function_value(ptId,i)*
00086 fN(fV.mapped_pts(ptId))*fV.JxW(ptId);
00087 }
00088
00089 }
