Hi, I'm Ataias Reis

Programmer, made in Brazil, Alma mater University of Brasilia

Code for Poisson with Dirichlet conditions using Explicit Method as well as python code to plot it

C++ code to get the solution

#include <iostream>
#include <chrono>
#include <functional>
#include <string>
#include <fstream>
#include <sstream>

using namespace std;

//compile: g++ -std=c++14 poissonDirichlet.cpp -o solver_dirichlet

// Thanks for an answer here, I am using this timer
//http://stackoverflow.com/questions/728068/how-to-calculate-a-time-difference-in-c
class Timer
{
public:
    Timer() : beg_(clock_::now()) {}
    void reset() { beg_ = clock_::now(); }
    double elapsed() const {
        return std::chrono::duration_cast<second_>
            (clock_::now() - beg_).count(); }

private:
    typedef std::chrono::high_resolution_clock clock_;
    typedef std::chrono::duration<double, std::ratio<1> > second_;
    std::chrono::time_point<clock_> beg_;
};


double abs(double x) {
  if(x < 0) return -x;
  return x;
}

// Remember we are considering number of points in x == y
double ** getMesh(int n) {
  // we use calloc to initialize everything with zeros
  double ** mesh = (double **) calloc (n, sizeof(double *));
  for (int i = 0; i < n; i++) {
    mesh[i] = (double *) calloc (n, sizeof(double));
  }

  return mesh;
}

/**
  g is a function of (x,y) that is applied on the boundaries
*/
void computeBoundaries(double ** mesh, int n, function<double(double, double)> g) {
  int i, j;
  double dx = 1 / (n - 1);
  i = 0;
  for (j = 0; j < n; j++) mesh[i][j] = g(i * dx, j * dx);
  i = n - 1;
  for (j = 0; j < n; j++) mesh[i][j] = g(i * dx, j * dx);
  j = 0;
  for (i = 0; i < n; i++) mesh[i][j] = g(i * dx, j * dx);
  j = n - 1;
  for (i = 0; i < n; i++) mesh[i][j] = g(i * dx, j * dx);

  return;
}

/**
 f is a function that expects mesh and i, j and computes the value of
 that point. It also takes the absolute change reference and adds to it.
 The return of this function is the absolute change
*/
double computeInternalPoints(double ** mesh, int n, function<double(double**, int, int, double*)> f) {
  double absoluteChange = 0.0;
  // only internal points, that's why there is the loop goes from 0 to n - 2
  for (int i = 1; i < n - 1; i++) {
    for (int j = 1; j < n - 1; j++) {
      mesh[i][j] = f(mesh, i, j, &absoluteChange);
    }
  }

  return absoluteChange;
}

void getSolutionDirichlet(double** mesh, int n, function<double(double, double)> g, double epsilon) {

  double absoluteChange = 0.0;
  double dx = 1.0 / (n - 1);

  // nabla f = u
  auto u = [] (double x, double y)  {return x*x + y*y;};

  // f_ij function
  auto f = [u, dx](auto mesh, int i, int j, auto absoluteChange) {
    double old_value = mesh[i][j];
    double c = mesh[i-1][j] + mesh[i+1][j] + mesh[i][j-1] + mesh[i][j+1];
    mesh[i][j] = (c - u(i * dx, j * dx)* dx*dx) / 4.0;
    * absoluteChange += abs(old_value - mesh[i][j]);
    return mesh[i][j];
  };

  // we just need to compute the boundaries once
  computeBoundaries(mesh, n, g);
  int k = 0;
  Timer tmr;
  do {
    absoluteChange = computeInternalPoints(mesh, n, f);
    k++;
  } while (absoluteChange > epsilon);

  cout << "Finished computing! Latest absolute change was: " << absoluteChange << endl;
  cout << "Total number of iterations on mesh: " << k << endl;
  double time_elapsed = tmr.elapsed();
  cout << "Total time spent: " << time_elapsed << " seconds" << endl;
  cout << "Total time per iteration: " << time_elapsed / k << " seconds" << endl;
  return;
}

void writeDataToFile(string filename, double ** mesh, int n){
  ofstream file(filename, ios::out|ios::binary);
  file.write((char*) &n, sizeof(int));
  for (int i = 0; i < n; i++) file.write((char*) mesh[i], n * sizeof(double));
  cout << "File size: " << file.tellp() << " bytes" << endl;
  file.close();
  cout << "File successfully saved at: " << filename << endl;
}

void freeArray(double ** mesh, int n) {
  for (int i = 0; i < n; i++) free(mesh[i]);
  free(mesh);
}


int main(int argc, char* argv[]) {
  if (argc < 3) {
    cout << "Usage: solver N outputFilename" << endl;
    return 1;
  }

  int n = 200;
  istringstream ss(argv[1]);
  if (!(ss >> n)){
    cerr << "Invalid number: " << argv[1] << '\n';
    return 1;
  }

  string filename(argv[2]);
  // Boundary condition
  auto g = [](double x, double y) {
    if (abs(x) < 0.0001) return (y - 0.5) * (y - 0.5);
    else if (abs(x) > 0.999) return (y - 0.5) * (y - 0.5);
    else if (abs(y) < 0.0001) return (x - 0.5) * (x-0.5);
    else if (abs(y) > 0.999) return (x - 0.5) * (x-0.5);
    return 0.0;
  };


  double epsilon = 0.0001;
  double ** mesh = getMesh(n);
  getSolutionDirichlet(mesh, n, g, epsilon);
  writeDataToFile(filename, mesh, n);
  // freeing memory is important
  freeArray(mesh, n);

  return 0;
}

Python Code to plot the solution

#!/usr/bin/env python3

import numpy as np
import struct
import os
import argparse
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt

def read_matrix(filename):
    f = open(filename, 'rb')
    n = struct.unpack('i', f.read(4))[0]
    print("Started reading file")
    mesh = np.zeros((n,n), dtype=np.double)
    for i in range(0, n):
        for j in range(0, n):
            mesh[i,j] = struct.unpack('d',f.read(8))[0]
    f.close()
    print("Finished reading file")
    return n, mesh

def plot(n, mesh):
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    X = np.linspace(0, 1, n)
    Y = np.linspace(0, 1, n)
    X, Y = np.meshgrid(X, Y)
    R = np.sqrt(X**2 + Y**2)
    Z = mesh
    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                           linewidth=0, antialiased=False)

    ax.zaxis.set_major_locator(LinearLocator(10))
    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

    fig.colorbar(surf, shrink=0.5, aspect=5)

    plt.show()

def main():
    parser = argparse.ArgumentParser(description='Read 2D array of doubles and plot it considering x in (0,1) and y in (0,1).')
    parser.add_argument('--filename', dest='filename', type=str,
                    help='the binary file that is going to be used', required=True)
    args = parser.parse_args()
    n, mesh = read_matrix(args.filename)

    plot(n, mesh)

if __name__ == '__main__':
    main()
Tags
comments powered by Disqus