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

Vitelotte: Removed dummy node stuff. (replaced by in solver logic)

parent 87e62fac
......@@ -37,7 +37,12 @@ public:
};
public:
inline ElementBuilderBase() {}
inline ElementBuilderBase() : m_size(0) {}
inline void begin(const Mesh& mesh) { m_size = mesh.nodesSize(); }
unsigned end(const Mesh& /*mesh*/) { return matrixSize(); }
unsigned matrixSize() const { return m_size; }
inline void setRhs(const Mesh& /*mesh*/, Matrix& rhs,
SolverError* /*error*/=0) {
......@@ -47,6 +52,9 @@ public:
MatrixType matrixType(const Mesh& /*mesh*/) const {
return MATRIX_SPD;
}
protected:
unsigned m_size;
};
......
......@@ -75,6 +75,7 @@ FemSolver<_Mesh, _ElementBuilder>::build()
unsigned nCoefficients = 0;
m_error.resetStatus();
m_solved = false;
m_elementBuilder.begin(*m_mesh);
for(FaceIterator elem = m_mesh->facesBegin();
elem != m_mesh->facesEnd(); ++elem)
{
......@@ -95,16 +96,19 @@ FemSolver<_Mesh, _ElementBuilder>::build()
}
assert(it == coefficients.end());
unsigned size = m_elementBuilder.end(*m_mesh);
// We use nodesSize instead of nNode because we index nodes with Node::idx()
m_stiffnessMatrix.resize(m_mesh->nodesSize(), m_mesh->nodesSize());
m_stiffnessMatrix.resize(size, size);
m_stiffnessMatrix.setFromTriplets(
coefficients.begin(), coefficients.end());
// std::cout << "S:\n" << Matrix(m_stiffnessMatrix) << "\n";
m_type = m_elementBuilder.matrixType(*m_mesh);
// std::cout << ((m_type == ElementBuilder::MATRIX_SPD)?
// "SPD\n": "Symetric\n");
m_b.resize(m_mesh->nodesSize(), m_mesh->nCoeffs());
m_b.resize(size, m_mesh->nCoeffs());
m_elementBuilder.setRhs(*m_mesh, m_b, &m_error);
if(m_error.status() == SolverError::STATUS_ERROR)
{
......@@ -321,9 +325,12 @@ FemSolver<_Mesh, _ElementBuilder>::solve()
for(unsigned i = 0; i < nUnknowns; ++i)
{
assert(!m_mesh->isConstraint(Node(m_perm[i])));
m_mesh->value(Node(m_perm[i])) = m_x.row(i).
template cast<typename Mesh::Scalar>();
if(m_perm[i] < m_mesh->nodesSize())
{
assert(!m_mesh->isConstraint(Node(m_perm[i])));
m_mesh->value(Node(m_perm[i])) = m_x.row(i).
template cast<typename Mesh::Scalar>();
}
}
m_solved = true;
......
......@@ -65,9 +65,14 @@ protected:
typedef Eigen::Matrix<Scalar, 3, 1> Vector3;
typedef Eigen::Matrix<Scalar, 6, 1> Vector6;
typedef typename Mesh::Halfedge Halfedge;
typedef std::map<unsigned, Halfedge> PGCMap;
public:
inline FVElementBuilder(Scalar sigma = Scalar(.5));
void begin(const Mesh& mesh);
unsigned nCoefficients(const Mesh& mesh, Face element,
SolverError* error=0) const;
......@@ -83,7 +88,9 @@ public:
}
private:
using Base::m_size;
Scalar m_sigma;
PGCMap m_pgcMap;
};
......
......@@ -22,6 +22,15 @@ FVElementBuilder<_Mesh, _Scalar>::FVElementBuilder(Scalar sigma)
{
}
template < class _Mesh, typename _Scalar >
void
FVElementBuilder<_Mesh, _Scalar>::begin(const Mesh& mesh) {
Base::begin(mesh);
m_pgcMap.clear();
}
template < class _Mesh, typename _Scalar >
unsigned
FVElementBuilder<_Mesh, _Scalar>::
......@@ -185,8 +194,12 @@ FVElementBuilder<_Mesh, _Scalar>::
0,
-elem.doubleArea()/(elem.edgeLength(1)*elem.edgeLength(2));
int ce1 = mesh.vertexGradientDummyNode(h1).idx();
int ce2 = mesh.vertexGradientDummyNode(h2).idx();
// int ce1 = mesh.vertexGradientDummyNode(h1).idx();
// int ce2 = mesh.vertexGradientDummyNode(h2).idx();
m_pgcMap.insert(std::make_pair(m_size, h1));
int ce1 = (m_size++);
m_pgcMap.insert(std::make_pair(m_size, h2));
int ce2 = (m_size++);
// std::cout << " ce1: " << ce1 << ", ce2: " << ce2 << "\n";
if(ce1 < 0 || ce2 < 0)
{
......@@ -213,24 +226,42 @@ FVElementBuilder<_Mesh, _Scalar>::
rhs.setZero();
for(typename Mesh::HalfedgeIterator hit = mesh.halfedgesBegin();
hit != mesh.halfedgesEnd(); ++hit) {
if(mesh.nVertexGradientConstraints(*hit) == 0)
continue;
typename Mesh::Vertex from = mesh.fromVertex(*hit);
typename Mesh::Vertex to = mesh. toVertex(*hit);
typename Mesh::Node n = mesh.vertexGradientDummyNode(*hit);
if(n.isValid()) {
bool v0c = mesh.isGradientConstraint(from);
const typename Mesh::Gradient& grad = mesh.gradientConstraint(v0c? from: to);
typename Mesh::Vector v = mesh.position(to) - mesh.position(from);
if(!v0c) v = -v;
typename Mesh::Value cons = grad * v;
rhs.row(n.idx()) = cons.template cast<Scalar>();
}
for(typename PGCMap::const_iterator it = m_pgcMap.begin();
it != m_pgcMap.end(); ++it)
{
unsigned index = it->first;
typename Mesh::Halfedge h = it->second;
typename Mesh::Vertex from = mesh.fromVertex(h);
typename Mesh::Vertex to = mesh. toVertex(h);
bool v0c = mesh.isGradientConstraint(from);
const typename Mesh::Gradient& grad = mesh.gradientConstraint(v0c? from: to);
typename Mesh::Vector v = mesh.position(to) - mesh.position(from);
if(!v0c) v = -v;
typename Mesh::Value cons = grad * v;
rhs.row(index) = cons.template cast<Scalar>();
}
// unsigned index = mesh.nodesCount();
// for(typename Mesh::HalfedgeIterator hit = mesh.halfedgesBegin();
// hit != mesh.halfedgesEnd(); ++hit) {
// if(mesh.nVertexGradientConstraints(*hit) == 0)
// continue;
// typename Mesh::Vertex from = mesh.fromVertex(*hit);
// typename Mesh::Vertex to = mesh. toVertex(*hit);
// typename Mesh::Node n = mesh.vertexGradientDummyNode(*hit);
// if(n.isValid()) {
// bool v0c = mesh.isGradientConstraint(from);
// const typename Mesh::Gradient& grad = mesh.gradientConstraint(v0c? from: to);
// typename Mesh::Vector v = mesh.position(to) - mesh.position(from);
// if(!v0c) v = -v;
// typename Mesh::Value cons = grad * v;
// rhs.row(n.idx()) = cons.template cast<Scalar>();
// }
// }
}
......
......@@ -441,10 +441,6 @@ public:
inline unsigned nVertexGradientConstraints(Halfedge h) const;
inline unsigned nVertexGradientConstraints(Face f) const;
inline Node vertexGradientDummyNode(Halfedge h) const
{ assert(hasVertexGradientConstraint());
return m_vertexGradientDummyNodes.at(h); }
/// \}
......@@ -678,8 +674,6 @@ protected:
HalfedgeProperty<Node> m_edgeValueNodes;
HalfedgeProperty<Node> m_edgeGradientNodes;
HalfedgeNodeMap m_vertexGradientDummyNodes;
unsigned m_deletedNodes;
NodeMatrix m_nodes;
NodeProperty<bool> m_ndeleted;
......
......@@ -178,14 +178,6 @@ VGMesh<_Scalar, _Dim, _Chan>::garbageCollection(unsigned flags)
m_vertexGradientConstraints.swap(vxGradConstraints);
// remap halfedges
HalfedgeNodeMap dummyNodes;
for(typename HalfedgeNodeMap::const_iterator hNode = m_vertexGradientDummyNodes.begin();
hNode != m_vertexGradientDummyNodes.end(); ++hNode)
{
dummyNodes.insert(std::make_pair(gcMap(hNode->first), hNode->second));
}
m_vertexGradientDummyNodes.swap(dummyNodes);
for(HalfedgeIterator hit = halfedgesBegin(); hit != halfedgesEnd(); ++hit)
{
for(unsigned ai = 0; ai < HALFEDGE_ATTRIB_COUNT; ++ai)
......@@ -388,23 +380,7 @@ void VGMesh<_Scalar, _Dim, _Chan>::
{
assert(hasVertexGradientConstraint());
bool newConstraint = !isGradientConstraint(v);
m_vertexGradientConstraints[v] = grad;
if(newConstraint)
{
HalfedgeAroundVertexCirculator hit = halfedges(v);
HalfedgeAroundVertexCirculator hend = hit;
do
{
if(!m_vertexGradientDummyNodes[*hit].isValid())
m_vertexGradientDummyNodes[*hit] = addNode();
if(!m_vertexGradientDummyNodes[oppositeHalfedge(*hit)].isValid())
m_vertexGradientDummyNodes[oppositeHalfedge(*hit)] = addNode();
++hit;
} while(hit != hend);
}
}
......@@ -414,15 +390,6 @@ void VGMesh<_Scalar, _Dim, _Chan>::removeGradientConstraint(Vertex v)
assert(hasVertexGradientConstraint());
m_vertexGradientConstraints.erase(v);
HalfedgeAroundVertexCirculator hit = halfedges(v);
HalfedgeAroundVertexCirculator hend = hit;
do
{
m_vertexGradientDummyNodes.erase(*hit);
m_vertexGradientDummyNodes.erase(oppositeHalfedge(*hit));
++hit;
} while(hit != hend);
}
......@@ -479,11 +446,6 @@ VGMesh<_Scalar, _Dim, _Chan>::deleteUnusedNodes()
m_ndeleted[edgeGradientNode(*hit)] = false;
}
}
for(typename HalfedgeNodeMap::const_iterator hNode = m_vertexGradientDummyNodes.begin();
hNode != m_vertexGradientDummyNodes.end(); ++hNode)
{
m_ndeleted[hNode->second] = false;
}
m_deletedNodes = 0;
for(unsigned ni = 0; ni < nodesSize(); ++ni)
{
......@@ -509,6 +471,7 @@ VGMesh<_Scalar, _Dim, _Chan>::addNode(const Eigen::DenseBase<Derived>& value)
return Node(nodesSize() - 1);
}
template < typename _Scalar, int _Dim, int _Chan >
bool
VGMesh<_Scalar, _Dim, _Chan>::hasUnknowns() const
......@@ -903,8 +866,6 @@ VGMesh<_Scalar, _Dim, _Chan>::copyVGMeshMembers(const Self& rhs)
m_edgeValueNodes = getHalfedgeProperty<Node>("h:edgeValueNode");
m_edgeGradientNodes = getHalfedgeProperty<Node>("h:edgeGradientNode");
m_vertexGradientDummyNodes = rhs.m_vertexGradientDummyNodes;
m_deletedNodes = rhs.m_deletedNodes;
m_nodes = rhs.m_nodes/*.template cast<Scalar>()*/;
m_ndeleted = nodeProperty<bool>("n:deleted");
......
......@@ -34,7 +34,6 @@ else:
exit(1)
out = open(filename, 'w')
sys.stdout = out
......@@ -69,23 +68,24 @@ for r in range(1, nrings+1):
outer += 1
if r == 1:
inner += 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")
print("pgc 0", .5/nrings, "0 0 0 0", .5/nrings, "0 0")
vx = 1 + 3 * nrings * (nrings - 1)
ring_size = nrings * 6
print("c ",
' '.join(map(lambda x: str(x + vx) + ' ' + str(x / ring_size), range(ring_size))),
' ', str(vx), ' 1', sep='')
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',
'0 1 0.5 0 1',
'0.166666 0.75 0.933 0 1',
'0.333333 0.25 0.933 0 1',
'.5 0 0.5 0 1',
'0.666666 0.25 0.067 0 1',
'0.833333 0.75 0.067 0 1',
'1 1 0.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='')
print("dcv 0", ' '.join(grad))
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