getmatrix.cpp

Go to the documentation of this file.

#include "getmatrix.h"
#include "userset.h"
 
#include <algorithm>
#include <cassert>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <list>
#include <vector>
using namespace std;
 
 
void Get_Matrix_Pattern(int const nelem, int const ndof_e, int const *const ia,
                        int &nnz, int* &id, int* &ik, double* &sk)
{
//  Determine the number of nodes
    int nnode = *max_element(ia, ia + ndof_e * nelem);
    ++nnode;                                 // node numberng: 0 ... nnode-1
    assert(*min_element(ia, ia + ndof_e * nelem) == 0); // numbering starts with 0 ?
 
//  Collect for each node those nodes it is connected to (multiple entries)
    vector< list<int> > cc(nnode);             //  cc[i] is the  list of nodes a node i is connected to
    for (int i = 0; i < nelem; ++i)
    {
        //auto &l2g = new_elem[i].Get_l2g();
        const int idx = ndof_e * i;
        for (int k = 0; k < ndof_e; ++k)
        {
            list<int> &cck = cc.at(ia[idx + k]);
            cck.insert( cck.end(), ia + idx, ia + idx + ndof_e );
        }
    }
//  Delete the multiple entries
    nnz = 0;
    for (auto &it : cc)
    {
        it.sort();
        it.unique();
        nnz += it.size();
        // cout << it.size() << " :: "; copy(it->begin(),it->end(), ostream_iterator<int,char>(cout,"  ")); cout << endl;
    }
 
// CSR data allocation
    id = new int [nnode + 1];                  // Allocate memory for CSR row pointer
    ik = new int [nnz];                        // Allocate memory for CSR column index vector
 
//  copy CSR data
    id[0] = 0;                                 // begin of first row
    for (size_t i = 0; i < cc.size(); ++i)
    {
        //cout << i << "   " << nid.at(i) << endl;;
        const list<int> &ci = cc.at(i);
        const auto nci = static_cast<int>(ci.size());
        id[i + 1] = id[i] + nci; // begin of next line
        copy(ci.begin(), ci.end(), ik + id[i] );
    }
 
    assert(nnz == id[nnode]);
    sk = new double [nnz];                      // Allocate memory for CSR column index vector
    return;
}
 
 
 
//  general routine for lin. triangular elements
 
void CalcElem(int const ial[3], double const xc[], double ske[3][3], double fe[3])
//void CalcElem(const int* __restrict__ ial, const double* __restrict__ xc, double* __restrict__ ske[3], double* __restrict__ fe)
{
    const int  i1  = 2 * ial[0],   i2 = 2 * ial[1],   i3 = 2 * ial[2];
    const double x13 = xc[i3 + 0] - xc[i1 + 0],  y13 = xc[i3 + 1] - xc[i1 + 1],
                 x21 = xc[i1 + 0] - xc[i2 + 0],  y21 = xc[i1 + 1] - xc[i2 + 1],
                 x32 = xc[i2 + 0] - xc[i3 + 0],  y32 = xc[i2 + 1] - xc[i3 + 1];
    const double jac = fabs(x21 * y13 - x13 * y21);
 
    ske[0][0] = 0.5 / jac * (y32 * y32 + x32 * x32);
    ske[0][1] = 0.5 / jac * (y13 * y32 + x13 * x32);
    ske[0][2] = 0.5 / jac * (y21 * y32 + x21 * x32);
    ske[1][0] = ske[0][1];
    ske[1][1] = 0.5 / jac * (y13 * y13 + x13 * x13);
    ske[1][2] = 0.5 / jac * (y21 * y13 + x21 * x13);
    ske[2][0] = ske[0][2];
    ske[2][1] = ske[1][2];
    ske[2][2] = 0.5 / jac * (y21 * y21 + x21 * x21);
 
    const double xm    = (xc[i1 + 0] + xc[i2 + 0] + xc[i3 + 0]) / 3.0,
                 ym    = (xc[i1 + 1] + xc[i2 + 1] + xc[i3 + 1]) / 3.0;
    fe[0] = fe[1] = fe[2] = 0.5 * jac * FunctF(xm, ym) / 3.0;
}
 
 
// general routine for lin. triangular elements,
// non-symm. matrix
// node numbering in element: lex. forward
 
void AddElem(int const ial[3], double const ske[3][3], double const fe[3],
             int const id[], int const ik[], double sk[], double f[])
{
    for (int i = 0; i < 3; ++i)
    {
        const int ii  = ial[i],           // row ii (global index)
                  id1 = id[ii],           // start and
                  id2 = id[ii + 1];       // end of row ii in matrix
        int ip  = id1;
        for (int j = 0; j < 3; ++j)       // no symmetry assumed
        {
            const int jj = ial[j];
            bool not_found = true;
            do       // find entry jj (global index) in row ii
            {
                not_found = (ik[ip] != jj);
                ++ip;
            }
            while (not_found && ip < id2);
 
#ifndef NDEBUG                 // compiler option -DNDEBUG switches off the check
            if (not_found)     // no entry found !!
            {
                cout << "Error in AddElem: (" << ii << "," << jj << ") ["
                     << ial[0] << "," << ial[1] << "," << ial[2] << "]\n";
                assert(!not_found);
            }
#endif
            sk[ip - 1] += ske[i][j];
        }
        f[ii] += fe[i];
    }
}
 
 
void DebugMatrix(int const nnode, int const id[], int const ik[], double const sk[])
{
//  ID points to first entry of row
//  no symmetry assumed
    cout << "\nMatrix  (nnz = " << id[nnode] << ")\n";
    int id1, id2 = -1;
    for (int row = 0; row < nnode; ++row)
    {
        cout << "Row " << row << " : ";
        id1 = id[row];
        id2 = id[row + 1];
        for (int j = id1; j < id2; ++j)
        {
            cout.setf(ios::right, ios::adjustfield);
            cout << "[" << setw(2) << ik[j] << "]  " << setw(4) << sk[j] << "  ";
        }
        cout << endl;
    }
    return;
}
 
 
void DebugVector(int const nnode, double const v[])
{
    cout << "\nVector  (nnode = " << nnode << ")\n";
    for (int j = 0; j < nnode; ++j)
    {
        cout.setf(ios::right, ios::adjustfield);
        cout << v[j] << "   ";
    }
    cout << endl;
 
    return;
}
 
// ----------------------------------------------------------------------------
 
 
int fetch(int const row, int const col, int const id[], int const ik[])
{
    int const id2 = id[row + 1];    // end   and
    int       ip  = id[row];        // start of recent row (global index)
 
    while (ip < id2 && ik[ip] != col)  // find index col (global index)
    {
        ++ip;
    }
#ifndef NDEBUG                 // compiler option -DNDEBUG switches off the check
    if (ip >= id2)
    {
        cout << "No column  " << col << "  in row  " << row << endl;
        ip = -1;
        assert(ip >= id2);
    }
#endif
    return ip;
}
 
//    w := f - K*u
void Defect(double w[], double const f[], double const u[],
            int const nnode, int const id[], int const ik[], double const sk[])
{
    for (int row = 0; row < nnode; ++row)
    {
        double wi = f[row];
        for (int ij = id[row]; ij < id[row + 1]; ++ij)
        {
            wi -= sk[ij] * u[ ik[ij] ];
        }
        w[row] = wi;
    }
    return;
}
 
void CrsMult(double w[], double const u[], int const nnode,
             int const id[], int const ik[], double const sk[])
{
    for (int row = 0; row < nnode; ++row)
    {
        double wi = 0.0;
        for (int ij = id[row]; ij < id[row + 1]; ++ij)
        {
            wi += sk[ij] * u[ ik[ij] ];
        }
        w[row] = wi;
    }
    return;
}
 
//----------------------------------------------------------------------
 
void GetDiag(int const nnode, int const id[], int const ik[], double const sk[], double d[])
{
    for (int row = 0; row < nnode; ++row)
    {
        const int ia = fetch(row, row, id, ik);
        assert(ia >= 0);
        d[row] = sk[ia];
    }
    return;
}
 
 
 
/*
        void GetMatrix(...)
 
        Computes matrix for Poisson eqution.
        We use a general code to calculate element matrices and
                local stiffness matrix K, no symmetry assumed
*/
void GetMatrix (int const nelem, int const ndof_e, int const ia[], int const nnode, double const xc[],
                int const nnz, int const id[], int const ik[], double sk[], double f[])
{
    assert(ndof_e == 3);
    assert(nnz == id[nnode]);
 
    for (int k = 0; k < id[nnode]; ++k)
    {
        sk[k] = 0.0;
    }
    for (int k = 0; k < nnode; ++k)
    {
        f[k] = 0.0;
    }
 
    double ske[3][3], fe[3];
    //  Loop over all elements
    for (int i = 0; i < nelem; ++i)
    {
        CalcElem(ia + 3 * i, xc, ske, fe);
        AddElem(ia + 3 * i, ske, fe, id, ik, sk, f);
    }
 
    //DebugMatrix(nnode, id, ik, sk);
 
    return;
}
 
/*  Apply b.c. defined in u on matrix and right hand side       */
/*  via penalty method.                                         */
// no symmetry of matrix
 
void ApplyDirichletBC(int const nx, int const ny, int const neigh[],
                      double const u[], int const id[], int const ik[] , double sk[], double f[])
{
    const double PENALTY = 1e6;
    const int dx = 1,
              dy = nx + 1,
              bl  = 0,
              br  = nx,
              tl  = ny * (nx + 1),
              tr  = (ny + 1) * (nx + 1) - 1;
    const int start[4] = { bl, br, tl, bl},
                         end[4] = { br, tr, tr, tl},
                                  step[4] = { dx, dy, dx, dy};
 
    for (int j = 0; j < 4; j++)
    {
        if (neigh[j] < 0)
        {
            for (int i = start[j]; i <= end[j]; i += step[j])
            {
                const int id1 = fetch(i, i, id, ik);
                assert(id1 >= 0);
//                if ( j!=1 )
                {
                    sk[id1] += PENALTY;         /* matrix weighted scaling feasible */
                    f[i] += PENALTY * u[i];
                }
//                else
//                {
//                   const double alfa = 1e-0;   // Robin B.C. with g=1e-2
//                   sk[id1] += alfa;
//                   f[i] += alfa*1e-2;
//                }
            }
        }
    }
    return;
}