Skip to content
Snippets Groups Projects
Commit 09891a78 authored by Oliver Sander's avatar Oliver Sander Committed by sander
Browse files

Remove some debugging code

[[Imported from SVN: r10106]]
parent 4c7ccc40
No related branches found
No related tags found
No related merge requests found
...@@ -277,32 +277,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve() ...@@ -277,32 +277,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
if (mgStep) if (mgStep)
corr = mgStep->getSol(); corr = mgStep->getSol();
//std::cout << "Correction: " << std::endl << corr_global << std::endl; //std::cout << "Correction: " << std::endl << corr << std::endl;
// Output correction for debugging
Dune::VTKWriter<typename GridType::LeafGridView> vtkWriter(grid_->leafGridView());
Dune::BlockVector<Dune::FieldVector<double,3> > displacement(x_.size());
for (size_t j=0; j<x_.size(); j++)
displacement[j] = x_[j] - identity_[j];
BasisType basis(grid_->leafGridView());
Dune::shared_ptr<VTKBasisGridFunction<BasisType,Dune::BlockVector<Dune::FieldVector<double,3> > > > vtkDisplacement
= Dune::make_shared<VTKBasisGridFunction<BasisType,Dune::BlockVector<Dune::FieldVector<double,3> > > >
(basis, displacement, "Displacement");
Dune::shared_ptr<VTKBasisGridFunction<BasisType,Dune::BlockVector<Dune::FieldVector<double,3> > > > vtkCorrection
= Dune::make_shared<VTKBasisGridFunction<BasisType,Dune::BlockVector<Dune::FieldVector<double,3> > > >
(basis, corr, "Correction");
Dune::shared_ptr<VTKBasisGridFunction<BasisType,Dune::BlockVector<Dune::FieldVector<double,3> > > > vtkGradient
= Dune::make_shared<VTKBasisGridFunction<BasisType,Dune::BlockVector<Dune::FieldVector<double,3> > > >
(basis, rhs, "Gradient");
vtkWriter.addVertexData(vtkDisplacement);
vtkWriter.addVertexData(vtkCorrection);
vtkWriter.addVertexData(vtkGradient);
vtkWriter.write("hencky_correction_" + std::to_string(i+1));
if (this->verbosity_ == NumProc::FULL) if (this->verbosity_ == NumProc::FULL)
std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl; std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
......
...@@ -123,8 +123,6 @@ protected: ...@@ -123,8 +123,6 @@ protected:
/** \brief An L2-norm, really. The H1SemiNorm class is badly named */ /** \brief An L2-norm, really. The H1SemiNorm class is badly named */
std::shared_ptr<H1SemiNorm<CorrectionType> > l2Norm_; std::shared_ptr<H1SemiNorm<CorrectionType> > l2Norm_;
public:
VectorType identity_;
}; };
#include "trustregionsolver.cc" #include "trustregionsolver.cc"
......
...@@ -138,7 +138,7 @@ int main (int argc, char *argv[]) try ...@@ -138,7 +138,7 @@ int main (int argc, char *argv[]) try
GridView gridView = grid->leafGridView(); GridView gridView = grid->leafGridView();
// typedef P1NodalBasis<GridView,double> FEBasis; // typedef P1NodalBasis<GridView,double> FEBasis;
typedef P2NodalBasis<GridView,double> FEBasis; typedef P1NodalBasis<GridView,double> FEBasis;
FEBasis feBasis(gridView); FEBasis feBasis(gridView);
// ///////////////////////////////////////// // /////////////////////////////////////////
...@@ -208,12 +208,6 @@ int main (int argc, char *argv[]) try ...@@ -208,12 +208,6 @@ int main (int argc, char *argv[]) try
for (size_t i=0; i<x.size(); i++) for (size_t i=0; i<x.size(); i++)
x[i] = v[i]; x[i] = v[i];
lambda = std::string("lambda x: (") + parameterSet.get<std::string>("identity") + std::string(")");
PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonIdentity(Python::evaluate(lambda));
SolutionType identity;
Functions::interpolate(feBasis, identity, pythonIdentity);
//////////////////////////////////////////////////////// ////////////////////////////////////////////////////////
// Main homotopy loop // Main homotopy loop
//////////////////////////////////////////////////////// ////////////////////////////////////////////////////////
...@@ -292,8 +286,6 @@ int main (int argc, char *argv[]) try ...@@ -292,8 +286,6 @@ int main (int argc, char *argv[]) try
baseTolerance baseTolerance
); );
solver.identity_ = identity;
//////////////////////////////////////////////////////// ////////////////////////////////////////////////////////
// Set Dirichlet values // Set Dirichlet values
//////////////////////////////////////////////////////// ////////////////////////////////////////////////////////
...@@ -333,6 +325,9 @@ int main (int argc, char *argv[]) try ...@@ -333,6 +325,9 @@ int main (int argc, char *argv[]) try
displacement[idx] = x[idx] - it->geometry().corner(0); displacement[idx] = x[idx] - it->geometry().corner(0);
} }
/////////////////////////////////
// Output result
/////////////////////////////////
Dune::shared_ptr<VTKBasisGridFunction<FEBasis,BlockVector<FieldVector<double,3> > > > vtkDisplacement Dune::shared_ptr<VTKBasisGridFunction<FEBasis,BlockVector<FieldVector<double,3> > > > vtkDisplacement
= Dune::make_shared<VTKBasisGridFunction<FEBasis,BlockVector<FieldVector<double,3> > > > = Dune::make_shared<VTKBasisGridFunction<FEBasis,BlockVector<FieldVector<double,3> > > >
...@@ -342,28 +337,6 @@ int main (int argc, char *argv[]) try ...@@ -342,28 +337,6 @@ int main (int argc, char *argv[]) try
} }
// //////////////////////////////
// Output result
// //////////////////////////////
// finally: compute the average deformation of the Neumann boundary
// That is what we need for the locking tests
FieldVector<double,3> averageDef(0);
for (size_t i=0; i<x.size(); i++)
if (neumannNodes[i][0])
{
averageDef += x[i];
std::cout << "i: " << i << ", pos: " << x[i] << std::endl;
}
averageDef /= neumannNodes.count();
if (mpiHelper.rank()==0)
{
std::cout << "Neumann value = " << parameterSet.get<std::string>("neumannValues") << std::endl;
std::cout << "Neumann value = " << parameterSet.get<FieldVector<double,dim> >("neumannValues") << " "
<< ", average deflection: " << averageDef << std::endl;
}
} catch (Exception e) { } catch (Exception e) {
std::cout << e << std::endl; std::cout << e << std::endl;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment