Skip to content
Snippets Groups Projects
Commit b6717460 authored by Oliver Sander's avatar Oliver Sander Committed by sander@PCPOOL.MI.FU-BERLIN.DE
Browse files

cleanup 

[[Imported from SVN: r9394]]
parent d0f18a18
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
test/checksymmetry.hh deleted 100644 → 0
+ 0
− 76
View file @ d0f18a18
#ifndef CHECK_SYMMETRY_HH
#define CHECK_SYMMETRY_HH
template <int blocksize>
void checkSymmetry(const Dune::BCRSMatrix<Dune::FieldMatrix<double,blocksize,blocksize> >& matrix,
const unsigned int* perm)
{
typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<double,blocksize,blocksize> >::row_type::ConstIterator ConstColumnIterator;
for (int i=0; i<matrix.N(); i++) {
ConstColumnIterator cIt = matrix[i].begin();
ConstColumnIterator cEndIt = matrix[i].end();
for (; cIt!=cEndIt; ++cIt) {
const Dune::FieldMatrix<double,blocksize,blocksize>& otherBlock = matrix[perm[i]][perm[cIt.index()]];
// Currently only for blocksize==2
Dune::FieldMatrix<double,blocksize,blocksize> diff, average;
diff[0][0] = (*cIt)[0][0] - otherBlock[0][0];
diff[0][1] = (*cIt)[0][1] + otherBlock[0][1];
diff[1][0] = (*cIt)[1][0] + otherBlock[1][0];
diff[1][1] = (*cIt)[1][1] - otherBlock[1][1];
average[0][0] = 0.5*((*cIt)[0][0] + otherBlock[0][0]);
average[0][1] = 0.5*((*cIt)[0][1] - otherBlock[0][1]);
average[1][0] = 0.5*((*cIt)[1][0] - otherBlock[1][0]);
average[1][1] = 0.5*((*cIt)[1][1] + otherBlock[1][1]);
if (diff.infinity_norm()/average.infinity_norm() > 1e-5) {
std::cout << "Asymmetry in matrix found!\n";
std::cout << "<<<< " << i << ", " << cIt.index() << std::endl << *cIt << std::endl;
std::cout << ">>>> " << perm[i] << ", " << perm[cIt.index()] << std::endl << otherBlock << std::endl;
exit(0);
}
}
}
}
template <int blocksize>
void checkSymmetry(const Dune::BlockVector<Dune::FieldVector<double,blocksize> >& vector,
const unsigned int* perm)
{
for (int i=0; i<vector.size(); i++) {
const Dune::FieldVector<double,blocksize>& otherBlock = vector[perm[i]];
// Currently only for blocksize==2
Dune::FieldVector<double,blocksize> diff, average;
diff[0] = vector[i][0] + otherBlock[0];
diff[1] = vector[i][1] - otherBlock[1];
average[0] = 0.5*(vector[i][0] - otherBlock[0]);
average[1] = 0.5*(vector[i][1] + otherBlock[1]);
if (diff.infinity_norm()/average.infinity_norm() > 1e-5) {
std::cout << "Asymmetry in vector found!\n";
std::cout << "<<<< " << i << std::endl << vector[i] << std::endl;
std::cout << ">>>> " << perm[i] << std::endl << otherBlock << std::endl;
exit(0);
}
}
}
#endif
test/finitedifferencecheck.hh deleted 100644 → 0
+ 0
− 87
View file @ d0f18a18
#ifndef FINITE_DIFFERENCE_CHECK_HH
#define FINITE_DIFFERENCE_CHECK_HH
template <class AssemblerType, class VectorType>
void finiteDifferenceCheck(AssemblerType& assembler, const VectorType& x)
{
VectorType rhs(x.size());
rhs = 0;
double eps = 1e-6;
const int blocksize = VectorType::block_type::size;
// ///////////////////////////////////////////////////////////
// Assemble gradient and hessian at the current position
// ///////////////////////////////////////////////////////////
assembler.assembleMatrix(x, rhs);
// ///////////////////////////////////////////////////////////
// Compare rhs with a finite-difference approximation
// ///////////////////////////////////////////////////////////
for (int i=0; i<rhs.size(); i++) {
for (int j=0; j<blocksize; j++) {
VectorType value0 = x;
VectorType value1 = x;
value0[i][j] -= eps;
value1[i][j] += eps;
double fd = (assembler.computeEnergy(value1) - assembler.computeEnergy(value0)) / (2*eps);
// relative error. NOTE: assembleMatrix currently assembles the NEGATIVE gradient!
double error = (-rhs[i][j] - fd)/(-rhs[i][j] + fd);
if (std::fabs(-rhs[i][j] - fd) > 1) {
std::cout << "ERROR! gradient: " << rhs[i][j] << ", fd: " << fd << std::endl;
exit(1);
}
//std::cout << "gradient: " << rhs[i][j] << ", fd: " << fd << std::endl;
}
}
//exit(0);
// ///////////////////////////////////////////////////////////
// Compare hessian with a finite-difference approximation
// ///////////////////////////////////////////////////////////
typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,blocksize,blocksize> > MatrixType;
const MatrixType& M = *assembler.getMatrix();
for (int i=0; i<M.N(); i++) {
typename MatrixType::row_type::ConstIterator cIt = M[i].begin();
typename MatrixType::row_type::ConstIterator cEndIt = M[i].end();
for (; cIt!=cEndIt; ++cIt) {
// Loop over all entries of the matrix block
for (int j=0; j<blocksize; j++) {
for (int k=0; k<blocksize; k++) {
// Is the matrix symmetric?
if ( std::abs((*cIt)[j][k] - M[cIt.index()][i][k][j]) > 1e-6) {
std::cout << "Matrix asymmetry: "
<< (*cIt)[j][k] << " -- " << M[cIt.index()][i][k][j] << std::endl;
//exit(1);
}
//assert( std::abs((*cIt)[j][k] - M[cIt.index()][i][k][j]) < 1e-6);
}
}
}
}
}
#endif
test/makesymmgrids.hh deleted 100644 → 0
+ 0
− 109
View file @ d0f18a18
#ifndef MAKE_SYMMETRIC_GRIDS_HH
#define MAKE_SYMMETRIC_GRIDS_HH
#include "../../common/arcofcircle.hh"
unsigned int perm_tire2d[18] = {16,
17,
14,
15,
12,
13,
10,
11,
8,
9,
6,
7,
4,
5,
2,
3,
0,
1};
void make2dTire(Dune::UGGrid<2,2>& grid, int numTangential, int numRadial)
{
using namespace Dune;
// Make tire
double innerRadius = 90;
double thickness = 10;
FieldVector<double,2> center;
center[0] = 0;
center[1] = 120;
grid.createbegin();
// ////////////////////////////////////////
// make parametrized boundary segments
// ////////////////////////////////////////
#if 1
std::vector<int> vertices(2);
for (int i=0; i<numTangential-1; i++) {
double lAngle = M_PI + M_PI*((double)i)/(numTangential-1);
double rAngle = M_PI + M_PI*((double)i+1)/(numTangential-1);
vertices[0] = i*numRadial; vertices[1] = (i+1)*numRadial;
grid.insertBoundarySegment(vertices, new ArcOfCircle(center, innerRadius, lAngle, rAngle));
vertices[0] = (i+1)*numRadial-1; vertices[1] = (i+2)*numRadial-1;
grid.insertBoundarySegment(vertices, new ArcOfCircle(center, innerRadius+thickness, lAngle, rAngle));
}
#endif
// ///////////////////////////
// make vertices
// ///////////////////////////
for (int i=0; i<numTangential; i++) {
double angle = M_PI*((double)i)/(numTangential-1);
for (int j=0; j<numRadial; j++) {
double radius = innerRadius + j*thickness/(numRadial-1);
FieldVector<double,2> v = center;
v[0] -= radius * std::cos(angle);
v[1] -= radius * std::sin(angle);
grid.insertVertex(v);
}
}
// ///////////////////////////
// make elements
// ///////////////////////////
std::vector<unsigned int> corners(4);
for (int i=0; i<numTangential-1; i++) {
for (int j=0; j<numRadial-1; j++) {
unsigned int base = i*numRadial + j;
corners[0] = base;
corners[1] = base+numRadial;
corners[2] = base+1;
corners[3] = base+numRadial+1;
grid.insertElement(cube, corners);
}
}
grid.createend();
}
#endif
test/neohooketest.cc deleted 100644 → 0
+ 0
− 52
View file @ d0f18a18
#include "config.h"
#include <dune/grid/uggrid.hh>
#include "../src/neohookeassembler.hh"
#include "makesymmgrids.hh"
#include "checksymmetry.hh"
#include "finitedifferencecheck.hh"
using namespace Dune;
int main() try
{
typedef BlockVector<FieldVector<double,2> > VectorType;
// Create the half-circle grid
typedef UGGrid<2,2> GridType;
GridType grid;
make2dTire(grid,9,2);
// Assemble
VectorType x(grid.size(0,2));
x = 0;
VectorType rhs(grid.size(0,2));
rhs = 0;
NeoHookeAssembler<GridType,1> assembler(grid);
assembler.assembleMatrix(x, rhs);
// Check matrix for symmetry
checkSymmetry(*assembler.getMatrix(), perm_tire2d);
// Check rhs for symmetry
checkSymmetry(rhs, perm_tire2d);
grid.globalRefine(2);
x.resize(grid.size(grid.maxLevel(),2));
x = 0;
for (int i=0; i<10; i++) {
x[i%x.size()][i%2] += i;
// Compare derivatives with a finite-difference approximation
finiteDifferenceCheck(assembler, x);
}
return 0;
} catch (Exception e) {
std::cout << e << std::endl;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment