1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#include "Deform2D.hpp"
#include <iostream>
#include <math.h>
#include <map>


extern "C" void openblas_set_num_threads(int);



namespace MeshKit {

Deform2D::Deform2D(vector<vector<double> > undeformed_cage_vertices, vector<vector<double> > deformed_cage_vertices)<--- Member variable 'Deform2D::num_interior' is not initialized in the constructor.
{
    num_vertices = undeformed_cage_vertices.size();
    //pass the data
	un_cage_vertices.insert(un_cage_vertices.begin(), undeformed_cage_vertices.begin(), undeformed_cage_vertices.end());
	cage_vertices.insert(cage_vertices.begin(), deformed_cage_vertices.begin(), deformed_cage_vertices.end());
}

void Deform2D::SetupInteriorNodes(vector<vector<double> > undeformed_in_nodes)<--- The function 'SetupInteriorNodes' is never used.
{
    num_interior = undeformed_in_nodes.size();
    
    un_InnerNodes.insert(un_InnerNodes.begin(), undeformed_in_nodes.begin(), undeformed_in_nodes.end());
	InnerNodes.insert(InnerNodes.begin(), undeformed_in_nodes.begin(), undeformed_in_nodes.end());    
}
void Deform2D::GetInteriorNodes(vector<vector<double> > &final_locations)<--- The function 'GetInteriorNodes' is never used.
{
   /*
   for (unsigned int i = 0; i < InnerNodes.size(); i++){
        for (unsigned int j = 0; j < InnerNodes[i].size(); j++){
            final_locations[i][j] = InnerNodes[i][j];
        }
   }
   */
   final_locations.clear();
   final_locations.insert(final_locations.begin(), InnerNodes.begin(), InnerNodes.end());
}

void Deform2D::Execute()
{
    const unsigned int vertex_size = un_cage_vertices.size();
    const unsigned int tmp_var = vertex_size + 3;<--- Variable 'tmp_var' is assigned a value that is never used.
    
    mat A = zeros<mat>(vertex_size+3,vertex_size+3);
    //mat::fixed<tmp_var,tmp_var> A;
    //A.zeros();
    /**************************/
    /*	    |	G	H |	 */
    /*	A = |		  |	 */			       
    /*	    |	H'	0 |      */			       
    /*****************************/
    //Assign H and H'
    for (unsigned int j = 0; j < vertex_size; j++){
		A(j,vertex_size) = 1.0;
		    A(j,vertex_size+1) = un_cage_vertices[j][0];
		A(j,vertex_size+2) = un_cage_vertices[j][1];
		A(vertex_size,j) = 1.0;
		A(vertex_size+1,j) = un_cage_vertices[j][0];
		A(vertex_size+2,j) = un_cage_vertices[j][1];
    }
	
    //assign G
    for (unsigned int i = 0; i < vertex_size; i++){
		for (unsigned int j = 0; j < vertex_size; j++){
			A(i,j) = TriHarmonicFun(un_cage_vertices[i], un_cage_vertices[j]);
		}
    }
    //done with the big matrix A
    //set up the displacement in x direction
    mat b=zeros<mat>(vertex_size+3,2);
    
    for (unsigned int i = 0; i < vertex_size; i++)
		b(i,0) = cage_vertices[i][0] - un_cage_vertices[i][0];
    //set up the displacement in y direction
    for (unsigned int i = 0; i < vertex_size; i++)
		b(i,1) = cage_vertices[i][1] - un_cage_vertices[i][1];

    std::cout << "Matrix size is " << vertex_size+3 << "*" << vertex_size+3 << "\n";
    openblas_set_num_threads(1);
    mat coeffs = solve(A,b);

    //compute the displacement in x and y direction for interior nodes
    for (unsigned int i = 0; i < num_interior; i++){
		InnerNodes[i][0] = un_InnerNodes[i][0] + InterpolatedFun(un_InnerNodes[i], coeffs.col(0));
		InnerNodes[i][1] = un_InnerNodes[i][1] + InterpolatedFun(un_InnerNodes[i], coeffs.col(1));
    }
    
}

double Deform2D::InterpolatedFun(vector<double> x, vec coeffs)
{
    double tmp = 0.0;
    double tmp_variable[3] = {1.0, x[0], x[1]};
    for (unsigned int i = 0; i < num_vertices; i++)
	    tmp += coeffs(i,0)*TriHarmonicFun(x, un_cage_vertices[i]);
    for (int i = 0; i < 3; i++)
	    tmp += coeffs(num_vertices+i,0)*tmp_variable[i];

    return tmp;
}

double Deform2D::TriHarmonicFun(vector<double> xi, vector<double> xj)
{
    double value = 0.0;
    double delta = sqrt(pow((xi[0] - xj[0]),2) + pow((xi[1] - xj[1]),2));
    value = delta*delta*log10(delta);<--- Variable 'value' is reassigned a value before the old one has been used.
    return value;
}

Deform2D::~Deform2D()
{
	cout << "It is over now in computing 2D deformation" << endl;
}

}