Commit c042d7c0 authored by Pit Nestle's avatar Pit Nestle

distortion_problem fixed, seems to work properly now.

parent db17e69e
Pipeline #7297 passed with stages
in 7 minutes and 56 seconds
......@@ -2,6 +2,7 @@
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
#include <iostream>
#include <cmath>
#if (COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
......@@ -18,30 +19,144 @@ class Problem_Distortion : public COMISO::NProblemInterfaceADOLC
{
public:
typedef Eigen::Matrix<adouble, 2, 2> Matrix2ad;
std::vector<Eigen::Vector3d> vertices;
Problem_Distortion() : COMISO::NProblemInterfaceADOLC()
typedef Eigen::Matrix<adouble, 2, 1> Vector2ad;
//Base Mesh
std::vector<Eigen::Vector3d> mesh_vertices;
int mesh_vertex_count;
//Parametrization Mesh
std::vector<Eigen::Vector2d> param_vertices;
int param_vertex_count;
//Triangles
struct Triangle { int a,b,c; };
std::vector<Triangle> triangles;
//Precalculated local Triangle Vectors and Area
std::vector<Eigen::Matrix2d> mesh_local_triangle_vectors;
std::vector<double> mesh_triangle_area;
Problem_Distortion() : COMISO::NProblemInterfaceADOLC(),
mesh_vertex_count(0),
param_vertex_count(0)
{}
virtual int n_unknowns()
{
return 0;
return param_vertex_count*2;
}
virtual void initial_x(double* _x)
{
if (this->param_vertex_count != this->mesh_vertex_count)
{
std::cerr << "Problem_Distorion::initial_x() : different vertex counts" << std::endl;
return;
}
for (int i = 0; i < this->param_vertex_count; i++)
{
_x[i*2] = this->param_vertices[i](0);
_x[i*2 + 1] = this->param_vertices[i](1);
std::cout << "Init \n" << mesh_vertices[i] << std::endl;
}
//local coordinate calculation for each triangle of the base mesh
for (Triangle t : triangles)
{
Eigen::Vector3d ab = mesh_vertices[t.b] - mesh_vertices[t.a]; // Vector 0 -> 1
Eigen::Vector3d ac = mesh_vertices[t.c] - mesh_vertices[t.a]; // Vector 0 -> 2
Eigen::Matrix2d local;
//Length of each vector
double ab_n = ab.norm();
double ac_n = ac.norm();
//Angle between ab and ac
double angle = std::acos((ab.dot(ac)) / (ab_n * ac_n));
//local ab vector
double ab_l_x = ab_n;
double ab_l_y = 0;
//local ac vector
double ac_l_x = std::cos(angle) * ac_n;
double ac_l_y = std::sin(angle) * ac_n;
//Triangle Area
double area = ab_l_x * ac_l_y * 0.5;
local << ab_l_x, ac_l_x, ab_l_y, ac_l_y; // x1, x2, y1, y2
this->mesh_local_triangle_vectors.push_back(local);
this->mesh_triangle_area.push_back(area);
}
}
virtual adouble eval_f_adouble(const adouble* _x)
{
adouble result = 0;
for (int i = 0; i < triangles.size(); i++)
{
Triangle t = triangles[i];
Vector2ad a;
Vector2ad b;
Vector2ad c;
a << _x[t.a*2],
_x[t.a*2 + 1];
b << _x[t.b*2],
_x[t.b*2 + 1];
c << _x[t.c*2],
_x[t.c*2 + 1];
Vector2ad ab = b - a;
Vector2ad ac = c - a;
Matrix2ad r;
r << ab, ac;
Matrix2ad q;
q(0, 0) = mesh_local_triangle_vectors[i](0, 0);
q(0, 1) = mesh_local_triangle_vectors[i](0, 1);
q(1, 0) = mesh_local_triangle_vectors[i](1, 0);
q(1, 1) = mesh_local_triangle_vectors[i](1, 1);
Matrix2ad j = q*r.inverse();
adouble j_norm = j.norm();
adouble j_inverse_norm = j.inverse().norm();
//Symmetric Dirichlet
adouble energy = j_norm*j_norm + j_inverse_norm*j_inverse_norm;
//
energy = (energy - 4) * mesh_triangle_area[i];
result += energy;
}
std::cout << "Energy : " << result.value() << std::endl;
return result;
}
virtual void store_result(const double* _x)
{
for (int i = 0; i < this->n_unknowns(); i++)
{
std::cout << _x[i] << std::endl;
}
}
......@@ -50,18 +165,54 @@ public:
virtual bool sparse_gradient() { return false;}
virtual bool sparse_hessian () { return true;}
void addVertex(double x, double y, double z)
void addMeshVertex(double x, double y, double z)
{
vertices.push_back(Eigen::Vector3d(x,y,z));
this->mesh_vertices.push_back(Eigen::Vector3d(x,y,z));
this->mesh_vertex_count++;
}
void addParamVertex(double x, double y)
{
this->param_vertices.push_back(Eigen::Vector2d(x,y));
this->param_vertex_count++;
}
bool addTriangle(int v0, int v1, int v2)
{
if (mesh_vertex_count < v0 || mesh_vertex_count < v1 || mesh_vertex_count < v2)
{
std::cerr << "Problem_Distorion::addTriangle() : Vertex index out of range." << std::endl;
return false;
}
this->triangles.push_back({v0, v1, v2});
return true;
}
bool addTriangle()
};
// Example main
int main(void)
{
Problem_Distortion pd;
pd.addMeshVertex(0, 0, 0);
pd.addMeshVertex(4, 0, 0);
pd.addMeshVertex(0, 3, 0);
pd.addMeshVertex(3, 4, 0);
pd.addParamVertex(1, 1);
pd.addParamVertex(2, 1);
pd.addParamVertex(1, 2);
pd.addParamVertex(2, 2);
pd.addTriangle(0, 1, 2);
pd.addTriangle(1, 2, 3);
COMISO::NewtonSolver ns;
ns.solve(&pd);
return 0;
}
......
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
#include <iostream>
#if (COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
#include <CoMISo/Utils/StopWatch.hh>
#include <vector>
#include <CoMISo/NSolver/LinearConstraint.hh>
#include <CoMISo/NSolver/NPDerivativeChecker.hh>
#include <CoMISo/NSolver/NProblemInterfaceADOLC.hh>
#include <CoMISo/NSolver/IPOPTSolver.hh>
#include <CoMISo/NSolver/NewtonSolver.hh>
class Problem_Distortion : public COMISO::NProblemInterfaceADOLC
{
public:
typedef Eigen::Matrix<adouble, 2, 2> Matrix2ad;
//Base Mesh
std::vector<Eigen::Vector3d> mesh_vertices;
unsigned int mesh_vertex_count;
//Parametrization Mesh
std::vector<Eigen::Vector2d> param_vertices;
unsigned int param_vertex_count;
//Triangle Vertex Indices
std::vector<unsigned int[3]> triangles;
Problem_Distortion() : COMISO::NProblemInterfaceADOLC(), vertex_count(0)
{}
virtual int n_unknowns()
{
return 0;
}
virtual void initial_x(double* _x)
{
}
virtual adouble eval_f_adouble(const adouble* _x)
{
adouble result = 0;
return result;
}
virtual void store_result(const double* _x)
{
}
virtual bool constant_gradient() const { return false; }
virtual bool constant_hessian() const { return false; }
virtual bool sparse_gradient() { return false;}
virtual bool sparse_hessian () { return true;}
void addMeshVertex(double x, double y, double z)
{
this->mesh_vertices.push_back(Eigen::Vector3d(x,y,z));
this->vertex_count++;
}
void addParamVertex(double x, double y, double z)
{
this->param_vertices.push_back(Eigen::Vector3d(x,y,z));
this->param_evertex_count++;
}
bool addTriangle(unsigned int v0, unsigned int v1, unsigned int v2)
{
if (vertex_count < v0 || vertex_count < v1 || vertex_count < c2)
{
std::cerr << "Problem_Distorion::addTriangle() : Vertex index out of range." << std::endl;
}
}
};
// Example main
int main(void)
{
return 0;
}
#else
int main(void)
{
std::cerr << "Warning: Example cannot be executed since either EIGEN3 or ADOLC is not available..." << std::endl;
return 0;
}
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment