Commit ddc12eab authored by Simon Boyé's avatar Simon Boyé

Vitelotte: exeriments with point gradient constraints.

parent 5f75d611
......@@ -137,6 +137,7 @@ FemSolver<_Mesh, _ElementBuilder>::solve()
StiffnessMatrix mat;
mat/*.template selfadjointView<Eigen::Lower>()*/ =
m_stiffnessMatrix.template selfadjointView<Eigen::Lower>().twistedBy(permInv);
// m_stiffnessMatrix.twistedBy(permInv);
// compute RHS
unsigned nUnknowns = m_ranges.back();
......@@ -159,23 +160,25 @@ FemSolver<_Mesh, _ElementBuilder>::solve()
m_x.setZero();
// Eigen::Matrix<float, 4, 2> g;
// g << 2, .5,
// 2, .5,
// 2, .5,
// g << 1., 0,
// 0, 1.,
// 0, 0,
// 0, 0;
// g /= 4;
//// g.setZero();
// typename Mesh::Vertex vert(5);
// typename Mesh::Vertex vert(0);
// typename Mesh::HalfedgeAroundVertexCirculator hit = m_mesh->halfedges(vert);
// typename Mesh::HalfedgeAroundVertexCirculator end = hit;
// do
// {
// Eigen::Vector2f e = (m_mesh->position(m_mesh->toVertex(*hit))
// -m_mesh->position(m_mesh->fromVertex(*hit))).normalized();
// - m_mesh->position(m_mesh->fromVertex(*hit))).normalized();
// Eigen::Vector4f v = g * e;
// std::cerr << "e: " << e.transpose() << " -> " << v.transpose() << "\n";
// m_x.row(m_mesh->edgeGradientNode(*hit).idx()) = v.cast<Scalar>();
// unsigned ni = m_mesh->edgeGradientNode(*hit).idx();
// std::cerr << "e " << ni << ": " << e.transpose() << " -> " << v.transpose() << "\n";
// m_x.row(ni) = v.cast<Scalar>();
// ++hit;
// } while(hit != end);
......@@ -199,6 +202,7 @@ FemSolver<_Mesh, _ElementBuilder>::solve()
L.startVec(j);
for(typename StiffnessMatrix::InnerIterator it(mat, start + j); it; ++it)
{
// if(it.index() >= nUnknowns) continue;
if(it.index() < j || it.index() >= nUnknowns) continue;
L.insertBackByOuterInnerUnordered(j, it.index() - start) = it.value();
}
......
......@@ -676,6 +676,12 @@ FVElementBuilder<_Mesh, _Scalar>::
{
typename Mesh::HalfedgeAroundFaceCirculator hit = mesh.halfedges(element);
// typename Mesh::HalfedgeAroundFaceCirculator hend = hit;
// do ++hit;
// while(mesh.toVertex(*hit).idx() != 0 && hit != hend);
// bool flat = mesh.toVertex(*hit).idx() == 0;
bool flat = false;
bool orient[3];
Vector p[3];
--hit;
......@@ -739,7 +745,7 @@ FVElementBuilder<_Mesh, _Scalar>::
- dx2[i] * dy2[j]
- dx2[j] * dy2[i]);
Scalar value = quadPointValue.sum() * (elem.doubleArea() / 3);
Scalar value = quadPointValue.sum() * (elem.doubleArea() / 6);
EIGEN_ASM_COMMENT("MYEND");
......@@ -755,9 +761,48 @@ FVElementBuilder<_Mesh, _Scalar>::
smNew(j, i) = value;
}
}
if(flat)
{
typedef Eigen::Matrix<Scalar, 9, 1> Vector9;
Vector9 fde1;
Vector9 fde2;
fde1 <<
-1.0L/2.0L*(2*elem.doubleArea()*elem.dldn(0, 1) + elem.doubleArea()*elem.dldn(1, 1) + 7*elem.edgeLength(1))/(elem.edgeLength(1)*elem.edgeLength(2)),
(1.0L/2.0L)*(elem.doubleArea()*elem.dldn(0, 0) + 2*elem.doubleArea()*elem.dldn(1, 0) - elem.edgeLength(0))/(elem.edgeLength(0)*elem.edgeLength(2)),
(1.0L/2.0L)*elem.doubleArea()*(elem.dldn(0, 0)*elem.edgeLength(1) - elem.dldn(1, 1)*elem.edgeLength(0) + 2*elem.dldn(2, 0)*elem.edgeLength(1) - 2*elem.dldn(2, 1)*elem.edgeLength(0))/(elem.edgeLength(0)*elem.edgeLength(1)*elem.edgeLength(2)),
-4/elem.edgeLength(2),
4/elem.edgeLength(2),
4/elem.edgeLength(2),
elem.doubleArea()/(elem.edgeLength(0)*elem.edgeLength(2)),
-elem.doubleArea()/(elem.edgeLength(1)*elem.edgeLength(2)),
0;
fde2 <<
-1.0L/2.0L*(2*elem.doubleArea()*elem.dldn(0, 2) + elem.doubleArea()*elem.dldn(2, 2) + 7*elem.edgeLength(2))/(elem.edgeLength(1)*elem.edgeLength(2)),
(1.0L/2.0L)*elem.doubleArea()*(elem.dldn(0, 0)*elem.edgeLength(2) + 2*elem.dldn(1, 0)*elem.edgeLength(2) - 2*elem.dldn(1, 2)*elem.edgeLength(0) - elem.dldn(2, 2)*elem.edgeLength(0))/(elem.edgeLength(0)*elem.edgeLength(1)*elem.edgeLength(2)),
(1.0L/2.0L)*(elem.doubleArea()*elem.dldn(0, 0) + 2*elem.doubleArea()*elem.dldn(2, 0) - elem.edgeLength(0))/(elem.edgeLength(0)*elem.edgeLength(1)),
-4/elem.edgeLength(1),
4/elem.edgeLength(1),
4/elem.edgeLength(1),
elem.doubleArea()/(elem.edgeLength(0)*elem.edgeLength(1)),
0,
-elem.doubleArea()/(elem.edgeLength(1)*elem.edgeLength(2));
smNew.row(7) += fde1;
// smNew.col(7) = fde1;
smNew.row(8) += fde2;
// smNew.col(8) = fde2;
std::cout << "Flat elem:\n";
std::cout << " p0: " << elem.point(0).transpose() << "\n";
std::cout << " p1: " << elem.point(1).transpose() << "\n";
std::cout << " p2: " << elem.point(2).transpose() << "\n";
std::cout << " Stiffness matrix:\n" << smNew << "\n";
}
}
Matrix9& sm = (method & 2)? smNew: smOld;
for(size_t i = 0; i < 9; ++i)
{
for(size_t j = 0; j < 9; ++j)
......
......@@ -60,4 +60,19 @@ for r in range(1, nrings+1):
# print("f", inner+v, outer+v, outer+v+1)
# print("f", inner+v, outer+v+1, inner+v+1)
# print("f", inner+r-1, outer+r-1, outer+r)
print("pc x o;", .5/nrings, "0 0 0 0", .5/nrings, "0 0; 0")
vx = 1 + 3 * nrings * (nrings - 1)
ring_size = nrings * 6
grad = [
'0:1,.5,0,1',
'0.166666:.75,.933,0,1',
'0.333333:.25,.933,0,1',
'.5:0,.5,0,1',
'0.666666:.25,.067,0,1',
'0.833333:.75,.067,0,1',
'1:1,.5,0,1',
]
print("dc o x; ", ' '.join(grad), "; ",
' '.join(map(lambda x: str(x + vx) + ':' + str(x / ring_size), range(ring_size))),
' ', str(vx), ':1', sep='')
......@@ -218,7 +218,10 @@ void Document::solve()
std::cout << "Solver error: " << m_fvSolver.errorString() << "\n";
}
exportPlot("plot.obj");
exportPlot("plot0.obj", 0);
exportPlot("plot1.obj", 1);
exportPlot("plot2.obj", 2);
exportPlot("plot3.obj", 3);
}
emit meshUpdated();
......@@ -352,7 +355,7 @@ void Document::loadMesh(const std::string& filename)
m_mesh.setAttributes(Mesh::FV_FLAGS);
if(m_mesh.nPointConstraints() || m_mesh.nCurves())
m_mesh.setNodesFromCurves();
m_mesh.setNodesFromCurves(0);
updateBoundingBox();
setSelection(MeshSelection());
......@@ -399,11 +402,10 @@ void Document::openSaveFinalMeshDialog()
}
void Document::exportPlot(const std::string& filename)
void Document::exportPlot(const std::string& filename, unsigned layer)
{
std::ofstream out(filename.c_str());
unsigned layer = 0;
unsigned nCount = 0;
std::vector<int> indexFromNode(m_solvedMesh.nNodes(), -1);
for(Mesh::FaceIterator fit = m_solvedMesh.facesBegin();
......
......@@ -113,7 +113,7 @@ public slots:
void openSaveSourceMeshDialog();
void openSaveFinalMeshDialog();
void exportPlot(const std::string& filename);
void exportPlot(const std::string& filename, unsigned layer=3);
signals:
void selectionChanged();
......
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