ImplementDirichletBC.cpp

#include <iostream>
#include "include/TriMatrix.h"
#include "include/Vector.h"

using namespace std;

// Implements the Dirichlet boundary condition to garantuee
// unique weak solution of boundary value problem
//
// $u(x_i)=g_D(x_i)$
//
// INPUTS:
// i      ... Index of node $x_i$
// g      ... $g_D(x_i)$
//
// OUTPUTS:
// matrix ... global stiffness matrix $\underline{K_h}$
// vector ... global load vector $\underline{f_h}$
void ImplementDirichletBC(long i, double g, TriMatrix& matrix, Vector& vector)
{
	for(long j=0;j<matrix.n;j++)
	{
		if(j!=i)
		{
			vector.elements[j]-=matrix.get(j,i)*g;
			matrix.set(j,i,0.);
			matrix.set(i,j,0.);
		}
		else
		{
			vector.elements[j]=matrix.get(i,i)*g;
		}
	}
}

ImplementRobinBC.cpp

#include <iostream>
#include "./include/TriMatrix.h"
#include "./include/Vector.h"

using namespace std;

// Implements the Robin boundary condition to garantuee
// unique weak solution of boundary value problem
//
// $u'(x_i)=\alpha(x_i)(g_R(x_i) - u(x_i))$
//
// INPUTS:
// i      ... Index of node $x_i$
// g      ... $g_R(x_i)$
// alpha  ... $\alpha(x_i)$
//
// OUTPUTS:
// matrix ... global stiffness matrix $\underline{K_h}$
// vector ... global load vector $\underline{f_h}$
void ImplementRobinBC(long i, double g, double alpha, TriMatrix& matrix, Vector& vector)
{
	long si=matrix.n-1;

	// Right boundary
	if(i==si)
	{
		matrix.set(si,si,matrix.get(si,si)+alpha);
		vector.elements[vector.n-1]+=alpha*g;
	}
	// Left boundary
	else if(i==0)
	{
		matrix.set(0,0,matrix.get(0,0)+alpha);
		vector.elements[0]+=alpha*g;
	}
	else
	{
		cerr << "No boundary point" << endl;
	}

}