Skip to content
Snippets Groups Projects
Commit c52b2bf5 authored by Elias Pipping's avatar Elias Pipping Committed by Elias Pipping
Browse files

Remove faulty gauss-seidel sample

parent 63e014d5
No related branches found
No related tags found
No related merge requests found
......@@ -5,15 +5,11 @@ check_PROGRAMS = \
test-gradient-method
bin_PROGRAMS = \
gauss-seidel-sample \
one-body-sample
one_body_sample_SOURCES = \
one-body-sample.cc
gauss_seidel_sample_SOURCES = \
gauss-seidel-sample.cc gaussseidel.hh
test_python_SOURCES = \
test-python.cc
......
/* -*- mode:c++; mode:semantic -*- */
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include "nicefunction.hh"
#include "mynonlinearity.hh"
#include "gaussseidel.hh"
int main() {
size_t const dimension = 3;
typedef Dune::FieldMatrix<double, dimension, dimension> LocalMatrix;
typedef Dune::FieldVector<double, dimension> LocalVector;
size_t const nodes = 10;
typedef Dune::FieldMatrix<LocalMatrix, nodes, nodes> GlobalMatrix;
typedef Dune::FieldVector<LocalVector, nodes> GlobalVector;
GlobalMatrix A;
// semi-random data
for (size_t i = 0; i < nodes; ++i)
for (size_t k = 0; k < nodes; ++k)
for (size_t j = 0; j < dimension; ++j)
for (size_t l = 0; l < dimension; ++l)
A[i][k][j][l] = ((i == k) && (j == l) ? 1 : 0);
GlobalVector b;
// semi-random data
for (size_t i = 0; i < nodes; ++i)
for (size_t j = 0; j < dimension; ++j)
b[i][j] = 1;
Dune::SampleFunction const f = Dune::SampleFunction(); // semi-random function
Dune::MyNonlinearity<dimension> const phi(f);
GlobalVector unu(LocalVector(0));
// semi-random data
for (size_t m = 0; m < nodes; ++m)
for (size_t j = 0; j < dimension; ++j)
unu[m][j] = m + j;
GaussSeidel(unu, A, b, phi);
std::cout << std::endl << std::endl;
GaussSeidel(unu, A, b, phi);
return 0;
}
/* -*- mode:c++; mode:semantic -*- */
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include "nicefunction.hh"
#include "mynonlinearity.hh"
#include "samplefunctional.hh"
// Just for debugging
template <int dim, int nodes>
double J_eval(
Dune::FieldMatrix<Dune::FieldMatrix<double, dim, dim>, nodes, nodes> const &
A,
Dune::FieldVector<Dune::FieldVector<double, dim>, nodes> const &b,
Dune::MyNonlinearity<dim> phi,
Dune::FieldVector<Dune::FieldVector<double, dim>, nodes> const &v) {
double ret = 0;
for (size_t i = 0; i < nodes; ++i) {
for (size_t k = 0; k < nodes; ++k) {
Dune::FieldMatrix<double, dim, dim> const localA = A[i][k];
Dune::FieldVector<double, dim> tmp;
localA.mv(v[k], tmp);
ret += 1 / 2 * (tmp * v[i]);
}
ret -= b[i] * v[i];
ret +=
phi(v[i]); // an example of phi^i(v[i]) (phi is independent of i here)
}
return ret;
}
template <int dim, int nodes>
void GaussSeidel(
Dune::FieldVector<Dune::FieldVector<double, dim>, nodes> &unu,
Dune::FieldMatrix<Dune::FieldMatrix<double, dim, dim>, nodes, nodes> const &
A,
Dune::FieldVector<Dune::FieldVector<double, dim>, nodes> const &b,
Dune::MyNonlinearity<dim> const &phi) {
for (size_t m = 0; m < nodes; ++m) {
typename Dune::SampleFunctional<dim>::SmallVector const localb = b[m];
typename Dune::SampleFunctional<dim>::SmallMatrix const localA = A[m][m];
Dune::SampleFunctional<dim> J(localA, localb, phi);
typename Dune::SampleFunctional<dim>::SmallVector correction;
for (size_t local_iteration = 0; local_iteration < 2;
++local_iteration) { // 2 is random here
Dune::minimise(J, unu[m], correction);
unu[m] += correction;
}
std::cout << J_eval(A, b, phi, unu) << std::endl; // debugging
}
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment