Skip to content

Commit

Permalink
removed updated lagrangian from frustrum DC example;
Browse files Browse the repository at this point in the history 
small changes in static solver
  • Loading branch information
@hverhelst
hverhelst committed Jun 23, 2021
1 parent bd69764 commit 117e7e4
Show file tree
Hide file tree
Showing 3 changed files with 36 additions and 25 deletions.
37 changes: 24 additions & 13 deletions examples/benchmark_FrustrumDC.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,10 +46,11 @@ int main (int argc, char** argv)
index_t impl = 1; // 1= analytical, 2= generalized, 3= spectral

int result = 0;
index_t maxit = 50;
index_t maxit = 20;
// Arc length method options
real_t dL = 2e-2; // General arc length
real_t tol = 1e-6;
real_t dL = -1; // General arc length
real_t tolU = 1e-3;
real_t tolF = 1e-3;

std::string wn("data.csv");

Expand Down Expand Up @@ -112,6 +113,7 @@ int main (int argc, char** argv)
*/
if (testCase == 0)
{
if (dL==-1) dL = 5e-2;
displ.setValue(-dL,3);
BCs.addCondition(boundary::north, condition_type::dirichlet, 0, 0, false, 0 ); // unknown 2 - z
BCs.addCondition(boundary::north, condition_type::dirichlet, 0, 0, false, 1 ); // unknown 2 - z
Expand All @@ -136,6 +138,7 @@ int main (int argc, char** argv)
}
else if (testCase == 1)
{
if (dL==-1) dL = 4e-2;
displ.setValue(-dL,3);

BCs.addCondition(boundary::north, condition_type::dirichlet, &displ, 0, false, 2 ); // unknown 1 - y
Expand Down Expand Up @@ -298,40 +301,44 @@ int main (int argc, char** argv)
return assembler->rhs();
};

gsSparseMatrix<> matrix;
gsVector<> vector;
assembler = new gsThinShellAssembler<3, real_t, true >(mp_def,dbasis,BCs,force,materialMatrix);
assembler->assemble();
gsSparseMatrix<> matrix = assembler->matrix();
gsVector<> vector = assembler->rhs();

gsStaticSolver<real_t> staticSolver(matrix,vector,Jacobian,Residual);
gsOptionList solverOptions = staticSolver.options();
solverOptions.setInt("Verbose",true);
solverOptions.setInt("MaxIterations",maxit);
solverOptions.setReal("Tolerance",tol);
solverOptions.setReal("ToleranceU",tolU);
solverOptions.setReal("ToleranceF",tolF);
staticSolver.setOptions(solverOptions);

real_t dL0 = dL;
int reset = 0;
gsMultiPatch<> mp_def0 = mp_def;
real_t indicator;

real_t D = 0;

for (index_t k=0; k<step; k++)
{
gsInfo<<"Load step "<< k<<"\n";

assembler = new gsThinShellAssembler<3, real_t, true >(mp_def,dbasis,BCs,force,materialMatrix);

stopwatch.restart();
stopwatch2.restart();
assembler->assemble();
displ.setValue(D - dL,3);
assembler->updateBCs(BCs);

time += stopwatch.stop();

matrix = assembler->matrix();
vector = assembler->rhs();
solVector = staticSolver.solveNonlinear();
totaltime += stopwatch2.stop();

if (!staticSolver.converged())
{
dL = dL/2;
displ.setValue(-dL,3);
displ.setValue(D -dL,3);
reset = 1;
mp_def = mp_def0;
gsInfo<<"Iterations did not converge\n";
Expand Down Expand Up @@ -382,11 +389,12 @@ int main (int argc, char** argv)
if (reset!=1)
{
dL = dL0;
displ.setValue(-dL,3);
displ.setValue(D - dL,3);
}
reset = 0;

mp_def0 = mp_def;
D -= dL;

gsInfo<<"--------------------------------------------------------------------------------------------------------------\n";
}
Expand All @@ -396,6 +404,9 @@ int main (int argc, char** argv)
gsInfo<<"Total ellapsed assembly time: \t\t"<<time<<" s\n";
gsInfo<<"Total ellapsed solution time (incl. assembly): \t"<<totaltime<<" s\n";

delete assembler;
delete materialMatrix;

return result;
}

Expand Down
2 changes: 1 addition & 1 deletion gsStaticSolver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,7 @@ class gsStaticSolver
const std::function < gsSparseMatrix<T> ( gsVector<T> const & ) > m_nonlinear;
const std::function < gsVector<T> ( gsVector<T> const & ) > m_residual;

gsVector<T> m_solVec;
gsVector<T> m_solVec, m_DeltaU, m_deltaU;

mutable index_t m_verbose;
bool m_NL;
Expand Down
22 changes: 11 additions & 11 deletions gsStaticSolver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,8 @@ gsVector<T> gsStaticSolver<T>::solveNonlinear()
if (residual==0) residual=1;
T residual0 = residual;
T residualOld = residual;
gsVector<T> DeltaU = gsVector<T>::Zero(m_solVec.rows());
gsVector<T> deltaU = gsVector<T>::Zero(m_solVec.rows());
m_DeltaU = gsVector<T>::Zero(m_solVec.rows());
m_deltaU = gsVector<T>::Zero(m_solVec.rows());
gsSparseMatrix<T> jacMat;


Expand All @@ -66,17 +66,17 @@ gsVector<T> gsStaticSolver<T>::solveNonlinear()

for (m_iterations = 0; m_iterations != m_maxIterations; ++m_iterations)
{
jacMat = m_nonlinear(m_solVec+DeltaU);
jacMat = m_nonlinear(m_solVec+m_DeltaU);
if (m_verbose==2)
{
gsInfo<<"Matrix: \n"<<jacMat.toDense()<<"\n";
gsInfo<<"Vector: \n"<<resVec<<"\n";
}
m_solver.compute(jacMat);
deltaU = m_solver.solve(resVec); // this is the UPDATE
DeltaU += m_relax * deltaU;
m_deltaU = m_solver.solve(resVec); // this is the UPDATE
m_DeltaU += m_relax * m_deltaU;

resVec = m_residual(m_solVec+DeltaU);
resVec = m_residual(m_solVec+m_DeltaU);
residual = resVec.norm();

if (m_verbose>0)
Expand All @@ -85,20 +85,20 @@ gsVector<T> gsStaticSolver<T>::solveNonlinear()
gsInfo<<std::setw(4)<<std::left<<m_iterations;
gsInfo<<std::setw(17)<<std::left<<residual;
gsInfo<<std::setw(17)<<std::left<<residual/residual0;
gsInfo<<std::setw(17)<<std::left<<m_relax * deltaU.norm();
gsInfo<<std::setw(17)<<std::left<<m_relax * deltaU.norm()/DeltaU.norm();
gsInfo<<std::setw(17)<<std::left<<m_relax * deltaU.norm()/m_solVec.norm();
gsInfo<<std::setw(17)<<std::left<<m_relax * m_deltaU.norm();
gsInfo<<std::setw(17)<<std::left<<m_relax * m_deltaU.norm()/m_DeltaU.norm();
gsInfo<<std::setw(17)<<std::left<<m_relax * m_deltaU.norm()/m_solVec.norm();
gsInfo<<std::setw(17)<<std::left<<math::log10(residualOld/residual0);
gsInfo<<std::setw(17)<<std::left<<math::log10(residual/residual0);
gsInfo<<"\n";
}

residualOld = residual;

if (m_relax * deltaU.norm()/m_solVec.norm() < m_toleranceU && residual/residual0 < m_toleranceF)
if (m_relax * m_deltaU.norm()/m_solVec.norm() < m_toleranceU && residual/residual0 < m_toleranceF)
{
m_converged = true;
m_solVec+=DeltaU;
m_solVec+=m_DeltaU;
break;
}
else if (m_iterations+1 == m_maxIterations)
Expand Down

0 comments on commit 117e7e4

Please sign in to comment.