libshell/example/StaticSolve.h at f8ec2da20745b2db6278f9e33a974f4e6d77e7d8 · evouga/libshell · GitHub

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
#ifndef STATICSOLVE_H
#define STATICSOLVE_H

#include <Eigen/Core>
#include <vector>

#include "libshell/MaterialModel.h"
#include "libshell/MeshConnectivity.h"
#include "libshell/ElasticShell.h"

template <class SFF>
void takeOneStep(const LibShell::MeshConnectivity &mesh,
    Eigen::MatrixXd &curPos,
    Eigen::VectorXd &curEdgeDOFs,
    const LibShell::MaterialModel<SFF> &mat,
    const LibShell::RestState &restState,
    double &reg)
{

    int nverts = (int)curPos.rows();
    int nedges = mesh.nEdges();
    int nedgedofs = SFF::numExtraDOFs;

    int freeDOFs = 3 * nverts + nedgedofs * nedges;

    while (true)
    {
        Eigen::VectorXd derivative;
        std::vector<Eigen::Triplet<double> > hessian;

        double energy = LibShell::ElasticShell<SFF>::elasticEnergy(mesh, curPos, curEdgeDOFs, mat, restState, &derivative, &hessian);

        Eigen::SparseMatrix<double> H(freeDOFs, freeDOFs);
        H.setFromTriplets(hessian.begin(), hessian.end());

        Eigen::VectorXd force = -derivative;
        Eigen::SparseMatrix<double> I(freeDOFs, freeDOFs);
        I.setIdentity();
        H += reg * I;

        Eigen::VectorXd maxvals(freeDOFs);
        maxvals.setZero();
        for (int k = 0; k < H.outerSize(); ++k)
        {
            for (Eigen::SparseMatrix<double>::InnerIterator it(H, k); it; ++it)
            {
                maxvals[it.row()] = std::max(maxvals[it.row()], std::fabs(it.value()));
            }
        }
        std::vector<Eigen::Triplet<double> > Dcoeffs;
        for (int i = 0; i < freeDOFs; i++)
        {
            double val = (maxvals[i] == 0.0 ? 1.0 : 1.0 / std::sqrt(maxvals[i]));
            Dcoeffs.push_back({ i,i, val });
        }
        Eigen::SparseMatrix<double> D(freeDOFs, freeDOFs);
        D.setFromTriplets(Dcoeffs.begin(), Dcoeffs.end());

        Eigen::SparseMatrix<double> DHDT = D * H * D.transpose();

        std::cout << "solving, original force residual: " << force.norm() << std::endl;
        Eigen::SimplicialLLT<Eigen::SparseMatrix<double> > solver(DHDT);
        if (solver.info() == Eigen::Success)
        {
            Eigen::VectorXd rhs = D * force;
            Eigen::VectorXd descentDir = D * solver.solve(rhs);
            std::cout << "solved" << std::endl;
            Eigen::MatrixXd newPos = curPos;
            for (int i = 0; i < nverts; i++)
            {
                newPos.row(i) += descentDir.segment<3>(3 * i);
            }
            Eigen::VectorXd newEdgeDofs = curEdgeDOFs + descentDir.segment(3 * nverts, nedgedofs * nedges);


            double newenergy = LibShell::ElasticShell<SFF>::elasticEnergy(mesh, newPos, newEdgeDofs, mat, restState, &derivative, NULL);
            force = -derivative;

            double forceResidual = force.norm();

            if (newenergy <= energy)
            {
                std::cout << "Old energy: " << energy << " new energy: " << newenergy << " force residual " << forceResidual << " pos change " << descentDir.segment(0, 3 * nverts).norm() << " theta change " << descentDir.segment(3 * nverts, nedgedofs*nedges).norm() << " lambda " << reg << std::endl;
                curPos = newPos;
                curEdgeDOFs = newEdgeDofs;
                reg /= 2.0;
                break;
            }
            else
            {
                std::cout << "Not a descent direction; old energy: " << energy << " new energy: " << newenergy << " lambda now: " << 2.0*reg << std::endl;
            }
        }
        else
        {
            std::cout << "Matrix not positive-definite, lambda now " << reg << std::endl;
        }

        reg *= 2.0;

    }
}

#endif