Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include "AAdapt_SPRSizeField.hpp"
00008 #include "AlbPUMI_FMDBMeshStruct.hpp"
00009 #include "Epetra_Import.h"
00010
00011 #include "pumi.h"
00012 #include "pumi_mesh.h"
00013 #include "apfPUMI.h"
00014 #include "apfSPR.h"
00015
00016 AAdapt::SPRSizeField::SPRSizeField(const Teuchos::RCP<AlbPUMI::AbstractPUMIDiscretization>& disc) :
00017 comm(disc->getComm()),
00018 mesh(disc->getFMDBMeshStruct()->apfMesh),
00019 esa(disc->getStateArrays().elemStateArrays),
00020 elemGIDws(disc->getElemGIDws()) {
00021 }
00022
00023 AAdapt::SPRSizeField::
00024 ~SPRSizeField() {
00025 }
00026
00027 void
00028 AAdapt::SPRSizeField::computeError() {
00029
00030 if ( sv_name.length() > 0 )
00031 computeErrorFromStateVariable();
00032 else
00033 computeErrorFromRecoveredGradients();
00034
00035 }
00036
00037
00038 void
00039 AAdapt::SPRSizeField::setParams(const Epetra_Vector* sol, const Epetra_Vector* ovlp_sol,
00040 double element_size, double err_bound,
00041 const std::string state_var_name) {
00042
00043 solution = sol;
00044 ovlp_solution = ovlp_sol;
00045 sv_name = state_var_name;
00046 rel_err = err_bound;
00047 std::vector<int> dims;
00048 esa[0][sv_name].dimensions(dims);
00049 num_qp = dims[1];
00050 cub_degree = getCubatureDegree(num_qp);
00051
00052 }
00053
00054 double AAdapt::SPRSizeField::getValue(ma::Entity* v) {
00055 return apf::getScalar(field,v,0);
00056 }
00057
00058 int AAdapt::SPRSizeField::getCubatureDegree(int num_qp) {
00059 switch(num_qp) {
00060 case 1:
00061 return 1;
00062 case 4:
00063 return 2;
00064 case 5:
00065 return 3;
00066 default:
00067 fprintf(stderr,"Invalid cubature degree");
00068 }
00069 }
00070
00071 void
00072 AAdapt::SPRSizeField::getFieldFromStateVariable(apf::Field* eps) {
00073 apf::Numbering* en = mesh->findNumbering("apf_element_numbering");
00074 apf::MeshIterator* it = mesh->begin(mesh->getDimension());
00075 apf::MeshEntity* e;
00076 while ((e = mesh->iterate(it))) {
00077 int elemID = apf::getNumber(en,e,0,0);
00078 int ws = elemGIDws[elemID].ws;
00079 int lid = elemGIDws[elemID].LID;
00080 for (int qp=0; qp < num_qp; qp++) {
00081 apf::Matrix3x3 value;
00082 for (int i=0; i<3; i++) {
00083 for (int j=0; j<3; j++) {
00084 value[i][j] = esa[ws][sv_name](lid,qp,i,j);
00085 }
00086 }
00087 apf::setMatrix(eps,e,qp,value);
00088 }
00089 }
00090 }
00091
00092 void
00093 AAdapt::SPRSizeField::computeErrorFromRecoveredGradients() {
00094
00095 apf::Field* f = mesh->findField("solution");
00096 apf::Field* sol_grad = apf::getGradIPField(f,"sol_grad",cub_degree);
00097 field = apf::getSPRSizeField(sol_grad,rel_err);
00098 apf::destroyField(sol_grad);
00099
00100 }
00101
00102
00103 void
00104 AAdapt::SPRSizeField::computeErrorFromStateVariable() {
00105
00106 apf::Field* eps = apf::createIPField(mesh,"eps",apf::MATRIX,cub_degree);
00107 getFieldFromStateVariable(eps);
00108 field = apf::getSPRSizeField(eps,rel_err);
00109 apf::destroyField(eps);
00110
00111 }
