diff --git a/src/model/westergaard.cpp b/src/model/westergaard.cpp index 9e87b07..ec2cc26 100644 --- a/src/model/westergaard.cpp +++ b/src/model/westergaard.cpp @@ -1,211 +1,215 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "fft_plan_manager.hh" #include "influence.hh" #include "loop.hh" #include "model.hh" #include "grid_view.hh" #include "static_types.hh" #include "westergaard.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ template Westergaard::Westergaard(Model* model) : IntegralOperator(model) { // Copy sizes std::array sizes; std::copy(model->getDiscretization().begin(), model->getDiscretization().end(), sizes.begin()); // Copy boundary sizes auto sizes_it = sizes.begin(); std::array boundary_hermitian_sizes; // Ignoring first dimension if model is volumetric if (trait::dimension > trait::boundary_dimension) ++sizes_it; std::copy(sizes_it, sizes.end(), boundary_hermitian_sizes.begin()); boundary_hermitian_sizes = GridHermitian::hermitianDimensions( boundary_hermitian_sizes); constexpr UInt nb_components = trait::components; buffer.setNbComponents(nb_components); buffer.resize(boundary_hermitian_sizes); influence.setNbComponents(nb_components * nb_components); influence.resize(boundary_hermitian_sizes); initInfluence(); } /* -------------------------------------------------------------------------- */ template template void Westergaard::fourierApply(Functor func, GridBase& in, GridBase& out) const try { Grid& i = dynamic_cast(in); Grid& full_out = dynamic_cast(out); FFTPlanManager::get().createPlan(i, buffer).forwardTransform(); /// Applying influence func(buffer, influence); /// Backward fourier transform on boundary only if constexpr (bdim == dim) FFTPlanManager::get().createPlan(full_out, buffer).backwardTransform(); else { auto view = make_view(full_out, 0); FFTPlanManager::get() .createPlan(view, buffer) .backwardTransform(); } } catch (const std::bad_cast& e) { TAMAAS_EXCEPTION("Neumann and dirichlet types do not match model type"); } /* -------------------------------------------------------------------------- */ template template void Westergaard::initFromFunctor(Functor func) { // Compute wavevectors for influence auto wavevectors = FFTransform::template computeFrequencies( influence.sizes()); // Get boundary physical size VectorProxy domain( this->model->getSystemSize()[(bdim == dim) ? 0 : 1]); // Normalize wavevectors wavevectors *= 2 * M_PI; wavevectors /= domain; // Apply functor Loop::stridedLoop(func, wavevectors, influence); // Set fundamental mode to zero MatrixProxy mat(influence(0)); mat = 0; } /* -------------------------------------------------------------------------- */ #define NEUMANN_BASIC(type) \ template <> \ void Westergaard::initInfluence() { \ auto E_star = model->getHertzModulus(); \ auto basic = [E_star] CUDA_LAMBDA(VectorProxy && q, \ - Complex & k) { k = 2. / q.l2norm(); }; \ + Complex & k) { \ + k = 2. / (E_star * q.l2norm()); \ + }; \ initFromFunctor(basic); \ } #define DIRICHLET_BASIC(type) \ template <> \ void Westergaard::initInfluence() { \ auto E_star = model->getHertzModulus(); \ auto basic = [E_star] CUDA_LAMBDA(VectorProxy && q, \ - Complex & k) { k = q.l2norm() / 2; }; \ + Complex & k) { \ + k = E_star * q.l2norm() / 2; \ + }; \ initFromFunctor(basic); \ } NEUMANN_BASIC(model_type::basic_1d); NEUMANN_BASIC(model_type::basic_2d); DIRICHLET_BASIC(model_type::basic_1d); DIRICHLET_BASIC(model_type::basic_2d); #undef NEUMANN_BASIC #undef DIRICHLET_BASIC /* -------------------------------------------------------------------------- */ template <> void Westergaard::initInfluence() { auto E = model->getYoungModulus(); auto nu = model->getPoissonRatio(); const Complex I(0, 1); auto surface = [E, nu, I] CUDA_LAMBDA(VectorProxy && q, MatrixProxy&& F) { Real q_norm = q.l2norm(); q /= q_norm; F(0, 0) = 2 * (1 + nu) * (1 - nu * q(0) * q(0)); F(1, 1) = 2 * (1 + nu) * (1 - nu * q(1) * q(1)); F(2, 2) = 2 * (1 - nu * nu); F(0, 1) = F(1, 0) = -q(0) * q(1) * 2 * nu * (1 + nu); F(0, 2) = I * q(0) * (1 + nu) * (1. - 2. * nu); F(1, 2) = I * q(1) * (1 + nu) * (1 - 2 * nu); F(2, 0) = -F(0, 2); F(2, 1) = -F(1, 2); F *= 1. / (E * q_norm); }; initFromFunctor(surface); } /* -------------------------------------------------------------------------- */ template void Westergaard::initInfluence() { TAMAAS_EXCEPTION("the requested operator has not been implemented"); } /* ------------------------------------------------------------------------ */ template void Westergaard::apply(GridBase& input, GridBase& output) const { auto apply = [](decltype(buffer)& buffer, const decltype(influence)& influence) { Loop::stridedLoop([] CUDA_LAMBDA(VectorProxy && i, MatrixProxy && F) { i = F * i; }, buffer, influence); }; fourierApply(apply, input, output); } /* -------------------------------------------------------------------------- */ /* Template instanciation */ /* -------------------------------------------------------------------------- */ template class Westergaard; template class Westergaard; template class Westergaard; template class Westergaard; template class Westergaard; template class Westergaard; template class Westergaard; template class Westergaard; template class Westergaard; template class Westergaard; template class Westergaard; template class Westergaard; /* -------------------------------------------------------------------------- */ __END_TAMAAS__ /* -------------------------------------------------------------------------- */ diff --git a/src/solvers/polonsky_keer_rey.cpp b/src/solvers/polonsky_keer_rey.cpp index 409f6fc..77a9216 100644 --- a/src/solvers/polonsky_keer_rey.cpp +++ b/src/solvers/polonsky_keer_rey.cpp @@ -1,293 +1,295 @@ /** * @file * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2016-2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "polonsky_keer_rey.hh" #include "elastic_functional.hh" #include "fft_plan_manager.hh" #include "loop.hh" #include "model_type.hh" #include #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ PolonskyKeerRey::PolonskyKeerRey(Model& model, const GridBase& surface, Real tolerance, type variable_type, type constraint_type) : ContactSolver(model, surface, tolerance), variable_type(variable_type), constraint_type(constraint_type) { if (model.getType() != model_type::basic_1d && model.getType() != model_type::basic_2d) { TAMAAS_EXCEPTION( "Model type is not compatible with PolonskyKeerRey solver"); } auto pressure_ = &model.getTraction(); switch (variable_type) { case gap: { + model.getBEEngine().registerDirichlet(); primal = this->_gap.get(); dual = pressure_; this->functional.addFunctionalTerm(new functional::ElasticFunctionalGap( model.getBEEngine(), this->surface), true); break; } case pressure: { + model.getBEEngine().registerNeumann(); primal = pressure_; dual = this->_gap.get(); this->functional.addFunctionalTerm( new functional::ElasticFunctionalPressure(model.getBEEngine(), this->surface), true); break; } } search_direction = allocateGrid(model.getType(), model.getDiscretization(), model.getTraction().getNbComponents()); projected_search_direction = allocateGrid(model.getType(), model.getDiscretization(), model.getTraction().getNbComponents()); } /* -------------------------------------------------------------------------- */ Real PolonskyKeerRey::solve(Real target) { Real G = 0, Gold = 1, error = 0, error_norm = 1; UInt n = 0; bool delta = false; *search_direction = 0; //*dual = 0; // Printing column headers std::cout << std::setw(5) << "Iter" << " " << std::setw(15) << "Cost_f" << " " << std::setw(15) << "Error" << '\n' << std::fixed; // Setting uniform value if constraint if (constraint_type == variable_type && std::abs(primal->sum()) <= 0) *primal = target; else if (constraint_type == variable_type) *primal *= target / primal->mean(); else if (constraint_type != variable_type) *primal = this->surface_stddev; do { // Computing functional gradient functional.computeGradF(*primal, *dual); Real dbar = meanOnUnsaturated(*dual); // Enforcing dual constraint via gradient if (constraint_type != variable_type) { *dual += 2 * target + dbar; } else { // Centering dual on primal > 0 *dual -= dbar; } // Computing gradient norm G = computeSquaredNorm(*dual); // Updating search direction (conjugate gradient) updateSearchDirection(delta * G / Gold); Gold = G; // Computing critical step Real tau = computeCriticalStep(target); // Update primal variable delta = updatePrimal(tau); // We should scale to constraint if (constraint_type == variable_type) *primal *= target / primal->mean(); error = computeError() / error_norm; Real cost_f = functional.computeF(*primal, *dual); printState(n, cost_f, error); } while (error > this->tolerance && n++ < this->max_iterations); // Final update of dual variable functional.computeGradF(*primal, *dual); Real dual_min = dual->min(); *dual -= dual_min; // Primal is pressure: need to make sure gap is positive if (variable_type == pressure) { model.getDisplacement() = *dual; model.getDisplacement() += this->surface; } // Primal is gap: need to make sure dual is positive (pressure + adhesion) else { model.getDisplacement() = *primal; model.getDisplacement() += this->surface; this->model.solveDirichlet(); model.getTraction() -= dual_min; } return error; } /* -------------------------------------------------------------------------- */ /** * Computes \f$ \frac{1}{\mathrm{card}(\{p > 0\})} \sum_{\{p > 0\}}{f_i} \f$ */ Real PolonskyKeerRey::meanOnUnsaturated(const GridBase& field) const { // Sum on unsaturated Real sum = Loop::reduce( [] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f : 0; }, *primal, field); // Number of unsaturated points UInt n_unsat = Loop::reduce( [] CUDA_LAMBDA(Real & p) -> UInt { return (p > 0); }, *primal); return sum / n_unsat; } /* -------------------------------------------------------------------------- */ /** * Computes \f$ \sum_{\{p > 0\}}{f_i^2} \f$ */ Real PolonskyKeerRey::computeSquaredNorm(const GridBase& field) const { return Loop::reduce( [] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f * f : 0; }, *primal, field); } /* -------------------------------------------------------------------------- */ /** * Computes \f$ \tau = \frac{ \sum_{\{p > 0\}}{q_i't_i} }{ \sum_{\{p > 0\}}{r_i' * t_i} } \f$ */ Real PolonskyKeerRey::computeCriticalStep(Real target) { switch (variable_type) { case gap: model.getBEEngine().solveDirichlet(*search_direction, *projected_search_direction); break; case pressure: model.getBEEngine().solveNeumann(*search_direction, *projected_search_direction); break; } Real rbar = meanOnUnsaturated(*projected_search_direction); if (constraint_type == variable_type) *projected_search_direction -= rbar; else { *projected_search_direction += 2 * target + rbar; } Real num = Loop::reduce( [] CUDA_LAMBDA(Real & p, Real & q, Real & t) { return (p > 0) ? q * t : 0; }, *primal, *dual, *search_direction); Real denum = Loop::reduce( [] CUDA_LAMBDA(Real & p, Real & r, Real & t) { return (p > 0) ? r * t : 0; }, *primal, *projected_search_direction, *search_direction); return num / denum; } /* -------------------------------------------------------------------------- */ /** * Update steps: * 1. \f$\mathbf{p} = \mathbf{p} - \tau \mathbf{t} \f$ * 2. Truncate all \f$p\f$ negative * 3. For all points in \f$I_\mathrm{na} = \{p = 0 \land q < 0 \}\f$ do \f$p_i = * p_i - \tau q_i\f$ */ bool PolonskyKeerRey::updatePrimal(Real step) { UInt na_num = Loop::reduce( [step] CUDA_LAMBDA(Real & p, Real & q, Real & t) -> UInt { p -= step * t; // Updating primal if (p < 0) p = 0; // Truncating neg values if (p == 0 && q < 0) { // If non-admissible state p -= step * q; return 1; } else return 0; }, *primal, *dual, *search_direction); return na_num == 0; } /* -------------------------------------------------------------------------- */ /** * Error is based on \f$ \sum{p_i q_i} \f$ */ Real PolonskyKeerRey::computeError() { /// Making sure dual respects constraint *dual -= dual->min(); Real error = primal->dot(*dual); Real norm = 1; if (variable_type == pressure) norm = std::abs(primal->sum() * this->surface_stddev); else norm = std::abs(dual->sum() * this->surface_stddev); norm *= primal->getNbPoints(); return std::abs(error) / norm; } /* -------------------------------------------------------------------------- */ /** * Do \f$\mathbf{t} = \mathbf{q}' + \delta \frac{R}{R_\mathrm{old}}\mathbf{t} * \f$ */ void PolonskyKeerRey::updateSearchDirection(Real factor) { Loop::loop( [factor] CUDA_LAMBDA(Real & p, Real & q, Real & t) { t = (p > 0) ? q + factor * t : 0; }, *primal, *dual, *search_direction); } /* -------------------------------------------------------------------------- */ void PolonskyKeerRey::addFunctionalTerm(functional::Functional* func) { functional.addFunctionalTerm(func, false); } __END_TAMAAS__ /* -------------------------------------------------------------------------- */