← Back to branch summary

~ubuntu-branches/ubuntu/maverick/dolfin/maverick

« back to all changes in this revision

Viewing changes to dolfin/ale/HarmonicSmoothing.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Johannes Ring
  • Date: 2008-09-16 08:41:20 UTC
  • Revision ID: james.westby@ubuntu.com-20080916084120-i8k3u6lhx3mw3py3
Tags: upstream-0.9.2
Import upstream version 0.9.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
dolfin/ale/HarmonicSmoothing.cpp
 
1
// Copyright (C) 2008 Anders Logg.
 
2
// Licensed under the GNU LGPL Version 2.1.
 
3
//
 
4
// First added:  2008-08-11
 
5
// Last changed: 2008-09-13
 
6
 
 
7
#include <dolfin/fem/Assembler.h>
 
8
#include <dolfin/log/log.h>
 
9
#include <dolfin/la/Matrix.h>
 
10
#include <dolfin/la/Vector.h>
 
11
#include <dolfin/la/solve.h>
 
12
#include <dolfin/mesh/MeshData.h>
 
13
#include "Poisson1D.h"
 
14
#include "Poisson2D.h"
 
15
#include "Poisson3D.h"
 
16
#include "HarmonicSmoothing.h"
 
17
 
 
18
using namespace dolfin;
 
19
 
 
20
//-----------------------------------------------------------------------------
 
21
void HarmonicSmoothing::move(Mesh& mesh, Mesh& new_boundary)
 
22
{
 
23
  // Choose form and function space
 
24
  FunctionSpace* V = 0;
 
25
  Form* form = 0;
 
26
  const uint D = mesh.topology().dim();
 
27
  const uint d = mesh.geometry().dim();
 
28
  switch (D)
 
29
  {
 
30
  case 1:
 
31
    V    = new Poisson1DFunctionSpace(mesh);
 
32
    form = new Poisson1DBilinearForm(*V, *V);
 
33
    break;
 
34
  case 2:
 
35
    V    = new Poisson2DFunctionSpace(mesh);
 
36
    form = new Poisson2DBilinearForm(*V, *V);
 
37
    break;
 
38
  case 3:
 
39
    V    = new Poisson3DFunctionSpace(mesh);
 
40
    form = new Poisson3DBilinearForm(*V, *V);
 
41
    break;
 
42
  default:
 
43
    error("Illegal mesh dimension %d for harmonic mesh smoothing.", D);
 
44
  };
 
45
 
 
46
  // Assemble matrix
 
47
  Matrix A;
 
48
  Assembler::assemble(A, *form);
 
49
 
 
50
  // Initialize vector
 
51
  const uint N = mesh.numVertices();
 
52
  Vector b(N);
 
53
 
 
54
  // Get array of dofs for boundary vertices
 
55
  const MeshFunction<uint>* vertex_map = new_boundary.data().mesh_function("vertex map");
 
56
  dolfin_assert(vertex_map);
 
57
  const uint num_dofs = vertex_map->size();
 
58
  const uint* dofs = vertex_map->values();
 
59
 
 
60
  // Modify matrix (insert 1 on diagonal)
 
61
  A.ident(num_dofs, dofs);
 
62
  A.apply();
 
63
 
 
64
  // Solve system for each dimension
 
65
  double* values = new double[num_dofs];
 
66
  double* new_coordinates = new double[d*N];
 
67
  Vector x;
 
68
  for (uint dim = 0; dim < d; dim++)
 
69
  {
 
70
    // Get boundary coordinates
 
71
    for (uint i = 0; i < new_boundary.numVertices(); i++)
 
72
      values[i] = new_boundary.geometry().x(i, dim);
 
73
    
 
74
    // Modify right-hand side
 
75
    b.set(values, num_dofs, dofs);
 
76
    b.apply();
 
77
    
 
78
    // Solve system
 
79
    solve(A, x, b, gmres, amg_hypre);
 
80
 
 
81
    // Get new coordinates
 
82
    x.get(new_coordinates + dim*N);
 
83
  }
 
84
 
 
85
  // Modify mesh coordinates
 
86
  MeshGeometry& geometry = mesh.geometry();
 
87
  for (uint dim = 0; dim < d; dim++)
 
88
    for (uint i = 0; i < N; i++)
 
89
      geometry.set(i, dim, new_coordinates[dim*N + i]);
 
90
 
 
91
  // Clean up
 
92
  delete V;
 
93
  delete form;
 
94
  delete [] values;
 
95
  delete [] new_coordinates;
 
96
}
 
97
//-----------------------------------------------------------------------------