Hi all,I have recently extended the penalized continuous contact brick to nonmatching meshes. However it doesn't seem to be very robust at all, especially when the global system contains other bricks that are based on Lagrange multipliers, so I am trying to find out if I have some bug somewhere or if it is just a theoretically expected behaviour.
In the attached files there is a quite minimal example demonstrating one case with good and one case with bad convergence. The only difference is the kind of the dirichlet conditions, penalized in the first case, with multipliers in the second one.
The examples can be run with ./test_contact input_simple and ./test_contact_fail input_simple respectively.It would be nice to have a second opinion before spending more time on debugging this issue, so that I can proceed later with the implementation of the Alart Curnier contact bricks for nonmatching meshes.
Best regards Kostas
% -*- octave -*- E1 = 210E+3; % [N/mm²] Young's Modulus for ... E2 = 210E+3; % [N/mm²] Young's Modulus for 34CrNiMo6 NU1 = 0.3; % [-] Poisson's ratio for ... NU2 = 0.3; % [-] Poisson's ratio for 34CrNiMo6 DX = 0.15; % [mm] resolution in x direction DY = 0.15; % [mm] resolution in y direction DZMIN1 = 0.2; % [mm] minimum element size in z direction DZMIN2 = 0.2; % [mm] minimum element size in z direction LZ1 = 0.2; % [mm] total height of the first body LZ2 = 0.2; % [mm] total height of the second body NZ_PER_LAYER1 = 1; % number of elements per mesh layer NZ_PER_LAYER2 = 1; % number of elements per mesh layer SURFACE1 = "plane_50.dat"; % filename for the 1st roughness topography SURFACE2 = "simple_2.dat"; % filename for the 2nd roughness topography DELTA = [0.05]; % [mm] relative approach of the bodies in contact PERIODIC_REPETITIONS_X = 1; PERIODIC_REPETITIONS_Y = 1;
plane_50.dat
Description: Netscape Proxy Auto Config
simple_2.dat
Description: Netscape Proxy Auto Config
#include "getfem/getfem_regular_meshes.h"
#include "getfem/getfem_models.h"
#include "getfem/getfem_model_solvers.h"
#include <getfem/getfem_Coulomb_friction.h>
#include "getfem/getfem_export.h"
using bgeot::short_type;
using bgeot::size_type;
using bgeot::scalar_type;
using bgeot::base_small_vector;
typedef getfem::model_real_plain_vector plain_vector;
typedef getfem::model_real_sparse_matrix sparse_matrix;
typedef gmm::row_matrix<gmm::rsvector<scalar_type> > rsr_matrix;
void read_matrix_from_file(const std::string &fname, std::vector< std::vector<double> > &mat)
{
mat.clear();
std::ifstream inpfile(fname.c_str());
std::string strline;
while (std::getline(inpfile,strline)) {
unsigned int first = strline.find_first_not_of(" ");
if (first == std::string::npos ||
strline[first] == '#')
continue; // skip empty or commented out line
std::stringstream linestream(strline);
std::string strvalue;
std::vector<double> vecline(0);
while (std::getline(linestream, strvalue, ',')) {
double dvalue = std::atof(strvalue.c_str());
vecline.push_back(dvalue);
}
mat.push_back(vecline);
}
}
// mesh regions
enum { XP = 0,
XM = 1,
YP = 2,
YM = 3,
ZP = 4,
ZM = 5
};
class block
{
private:
int id;
double dx,dy,dz,lz;
double lambda,mu;
int nx,ny,nz,zsteps;
bool reversed;
double delta;
std::vector< std::vector<double> > roughnessmat;
std::vector<getfem::mesh> meshes;
std::vector<getfem::mesh_fem> mfus, mfds;
std::vector<getfem::mesh_im> mims;
public:
block(int id_, double dx_, double dy_, double dz_, double lz_, int nz_,
double E, double nu, const std::string &surface_file, bool uniform_dz=false)
: id(id_), dx(dx_), dy(dy_), dz(dz_), lz(lz_), nx(0), ny(0), nz(nz_),
reversed(false), delta(0),
meshes(0), mfus(0), mfds(0), mims(0)
{
zsteps = int(floor(log2(1+lz/(dz*nz))));
if (uniform_dz)
dz = lz/(nz*(pow(2,zsteps)-1));
lambda = E*nu / ((1+nu)*(1-2*nu));
mu = E/(2*(1+nu));
read_matrix_from_file(surface_file, roughnessmat);
meshes.resize(zsteps);
mfus.resize(zsteps);
mfds.resize(zsteps);
mims.resize(zsteps);
}
int nx_max(void) { return ny_max() >= 0 ? roughnessmat[0].size()-1 : 0; }
int ny_max(void) { return roughnessmat.size()-1; }
void build_meshes(int nx_, int ny_, int dir, int nnx=0, int nny=0)
{
nx = std::min(nx_, nx_max());
ny = std::min(ny_, ny_max());
if (nnx < 1 || nnx > nx)
nnx = 1;
if (nny < 1 || nny > ny)
nny = 1;
int nx0 = nx/nnx;
nx = nx0*nnx;
int ny0 = ny/nny;
ny = ny0*nny;
reversed = (dir < 0);
bgeot::base_node org;
std::vector<bgeot::base_small_vector> vect(3);
std::vector<int> ref(3);
for (int i=0; i < zsteps; i++) {
size_type nX = nnx * std::max(1,(int)round(nx0/pow(2.,i)));
size_type nY = nny * std::max(1,(int)round(ny0/pow(2.,i)));
ref[0] = nX; ref[1] = nY; ref[2] = nz;
double dX = (dx*nx)/nX;
double dY = (dy*ny)/nY;
double dZ = dz*pow(2,i);
double stretch_factor=1.;
if (i == zsteps-1) // for non-uniform dz
stretch_factor = (lz - (pow(2,i)-1)*dz*nz)/(dZ*nz);
vect[0] = bgeot::base_small_vector(dX, 0.0, 0.0);
vect[1] = bgeot::base_small_vector(0.0, dY, 0.0);
vect[2] = bgeot::base_small_vector(0.0, 0.0, stretch_factor*dZ);
if (dir >= 0)
org = bgeot::base_node(0.0, 0.0, (dZ-dz)*nz);
else
org = bgeot::base_node(0.0, 0.0, -((1+stretch_factor)*dZ-dz)*nz);
getfem::parallelepiped_regular_mesh(meshes[i], 3, org, vect.begin(), ref.begin());
// define regions on the boundary
getfem::mesh &m = meshes[i];
for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
int nb_faces = m.structure_of_convex(cv)->nb_faces();
for (short_type fc=0; fc < nb_faces; fc++) {
if (!m.is_convex_having_neighbour(cv,fc)) {
base_small_vector un = m.normal_of_face_of_convex(cv,fc,0);
gmm::scale(un, 1/gmm::vect_norm2(un));
if (gmm::abs(un[0] - 1.0) < 1.0E-3)
m.region(XP).add(cv,fc);
else if (gmm::abs(un[0] + 1.0) < 1.0E-3)
m.region(XM).add(cv,fc);
else if (gmm::abs(un[1] - 1.0) < 1.0E-3)
m.region(YP).add(cv,fc);
else if (gmm::abs(un[1] + 1.0) < 1.0E-3)
m.region(YM).add(cv,fc);
else if (gmm::abs(un[2] - 1.0) < 1.0E-3)
m.region(ZP).add(cv,fc);
else if (gmm::abs(un[2] + 1.0) < 1.0E-3)
m.region(ZM).add(cv,fc);
}
} // fc
} // cv
}
// apply surface roughness topography
for (dal::bv_visitor ip(meshes[0].points_index()); !ip.finished(); ++ip) {
for (int i=nz-1; i >= 0; i--) {
double z = ((dir >= 0) ? dz : -dz) * i;
double scale = ((dir >= 0) ? 1 : -1) * double(nz-i)/double(nz);
if (gmm::abs(meshes[0].points()[ip][2]-z) < 1e-10) {
size_type ix=round(meshes[0].points()[ip][0]/dx);
size_type iy=round(meshes[0].points()[ip][1]/dy);
meshes[0].points()[ip][2] -= scale * roughnessmat[iy][ix];
}
}
}
}
void define_mesh_fems(void)
{
// mesh_fem's and integration methods
for (int i=0; i < zsteps; i++) {
mfus[i].init_with_mesh(meshes[i]);
mfds[i].init_with_mesh(meshes[i]);
mfus[i].set_qdim(3);
mfds[i].set_qdim(1);
mfus[i].set_classical_finite_element(1);
mfds[i].set_classical_finite_element(1);
mims[i].init_with_mesh(meshes[i]);
mims[i].set_integration_method(getfem::int_method_descriptor("IM_GAUSS_PARALLELEPIPED(3,5)"));
}
}
void add_linear_elasticity_to_model(getfem::model &md)
{
md.add_initialized_scalar_data(lambda_str(), lambda);
md.add_initialized_scalar_data(mu_str(), mu);
for (int i=0; i < zsteps; i++) {
md.add_fem_variable(var_str(i), mfus[i]);
getfem::add_isotropic_linearized_elasticity_brick(md, mims[i], var_str(i), lambda_str(), mu_str());
}
}
void add_dirichlet_conditions_to_model(getfem::model &md, double disp)
{
delta = disp;
int i=zsteps-1;
int reg=(reversed ? ZM : ZP);
plain_vector W(mfus[i].nb_dof());
gmm::clear(W);
dal::bit_vector cn=mfus[i].dof_on_region(reg);
for (dal::bv_visitor dof(cn); !dof.finished(); ++dof)
if (dof % 3 == 2)
W[dof] = disp;
md.add_initialized_fem_data(dirichlet_displacements_str(), mfus[i], W);
// getfem::add_Dirichlet_condition_with_multipliers
// (md, mims[i], var_str(i), mfus[i], reg, dirichlet_displacements_str());
getfem::add_Dirichlet_condition_with_penalization
(md, mims[i], var_str(i), 1e9, reg, dirichlet_displacements_str());
}
void set_dirichlet_conditions_constant(getfem::model &md, double disp)
{
delta = disp;
int i=zsteps-1;
int reg=(reversed ? ZM : ZP);
plain_vector W(mfus[i].nb_dof());
gmm::clear(W);
dal::bit_vector cn=mfus[i].dof_on_region(reg);
for (dal::bv_visitor dof(cn); !dof.finished(); ++dof)
if (dof % 3 == 2)
W[dof] = disp;
gmm::copy(W, md.set_real_variable(dirichlet_displacements_str()));
}
void export_results_from_model(getfem::model &md, const std::string &suffix="")
{
for (int i=0; i < zsteps; i++) {
std::stringstream ss_id_i;
ss_id_i << id << "_" << i;
std::stringstream ss_fname;
ss_fname << "./output/roughcontact" << id << "_" << i << suffix << ".vtk";
getfem::vtk_export exp(ss_fname.str(), true);
exp.exporting(meshes[i]);
exp.write_point_data(mfus[i], md.real_variable(var_str(i)), "elastostatic_displacement" + ss_id_i.str());
plain_vector VonMisesStresses(mfds[i].nb_dof());
getfem::compute_isotropic_linearized_Von_Mises_or_Tresca
(md, var_str(i), lambda_str(), mu_str(), mfds[i], VonMisesStresses, false);
exp.write_point_data(mfds[i], VonMisesStresses, "von mises stresses");
}
}
const std::string lambda_str(void)
{
std::stringstream ss_lambda;
ss_lambda << "lambda" << id;
return ss_lambda.str();
}
const std::string mu_str(void)
{
std::stringstream ss_mu;
ss_mu << "mu" << id;
return ss_mu.str();
}
const std::string var_str(int i)
{
if (i >= 0 && i <= zsteps-1) {
std::stringstream ss_var;
ss_var << "u" << id << "_" << i;
return ss_var.str();
}
else
return "";
}
const std::string mortar_mult_str(int i)
{
if (i >= 1 && i <= zsteps-1) {
std::stringstream ss_mult;
ss_mult << "mortar_mult" << id << "_" << i;
return ss_mult.str();
}
else
return "";
}
const std::string periodicity_mult_str(void)
{
std::stringstream ss_mult;
ss_mult << "periodicity_mult" << id;
return ss_mult.str();
}
const std::string dirichlet_displacements_str()
{
std::stringstream ss_dirichlet(std::stringstream::out);
ss_dirichlet << "dirichlet_displacements" << id;
return ss_dirichlet.str();
}
const getfem::mesh_im& mim(int i)
{
return mims[i];
}
double suggested_augmentation_parameter(void)
{
return mu * (3*lambda + 2*mu) / (lambda + mu);
}
};
int main(int argc, char* argv[])
{
bgeot::md_param PARAM;
PARAM.read_command_line(argc, argv);
double dx = PARAM.real_value("DX","resolution in x direction [mm]");
double dy = PARAM.real_value("DY","resolution in y direction [mm]");
double lz1 = PARAM.real_value("LZ1","total height of the first body [mm]");
double lz2 = PARAM.real_value("LZ1","total height of the second body [mm]");
double dz1 = PARAM.real_value("DZMIN1","minimum element size in z direction (1) [mm]");
double dz2 = PARAM.real_value("DZMIN2","minimum element size in z direction (2) [mm]");
const std::vector<bgeot::md_param::param_value> &delta_params = PARAM.array_value("DELTA");
std::vector<double> deltas(delta_params.size());
for (size_type i=0; i < deltas.size(); ++i) {
GMM_ASSERT1(delta_params[i].type_of_param() == bgeot::md_param::REAL_VALUE,
"DELTA should be a scalar array.");
deltas[i] = delta_params[i].real();
}
int nz1 = PARAM.int_value("NZ_PER_LAYER1","number of elements per mesh layer (1)");
int nz2 = PARAM.int_value("NZ_PER_LAYER2","number of elements per mesh layer (2)");
double E1 = PARAM.real_value("E1","Young's Modulus (1) [N/mm²]");
double E2 = PARAM.real_value("E2","Young's Modulus (2) [N/mm²]");
double nu1 = PARAM.real_value("NU1","Poisson's ratio (1) [-]");
double nu2 = PARAM.real_value("NU2","Poisson's ratio (2) [-]");
std::string roughness1_file = PARAM.string_value("SURFACE1","filename for the 1st roughness topography");
std::string roughness2_file = PARAM.string_value("SURFACE2","filename for the 1st roughness topography");
block block1(1,dx,dy,dz1,lz1,nz1,E1,nu1,roughness1_file);
block block2(2,dx,dy,dz2,lz2,nz2,E2,nu2,roughness2_file);
size_type nx = std::min(block1.nx_max(), block2.nx_max());
size_type ny = std::min(block1.ny_max(), block2.ny_max());
int nnx = PARAM.int_value("PERIODIC_REPETITIONS_X");
int nny = PARAM.int_value("PERIODIC_REPETITIONS_Y");
// meshes
block1.build_meshes(nx,ny,1,nnx,nny);
block2.build_meshes(nx,ny,-1,nnx,nny);
// mesh_fem's
block1.define_mesh_fems();
block2.define_mesh_fems();
// model
getfem::model md;
block1.add_linear_elasticity_to_model(md);
block2.add_linear_elasticity_to_model(md);
block1.add_dirichlet_conditions_to_model(md,-deltas[0]/2);
block2.add_dirichlet_conditions_to_model(md,deltas[0]/2);
md.add_initialized_scalar_data("r", 1e6);
getfem::add_penalized_contact_between_nonmatching_meshes_brick
(md, block1.mim(0), block1.var_str(0), block2.var_str(0), "r", ZM, ZP);
std::cout << "Solving model with " << md.nb_dof() << " degrees of freedom at "
<< deltas.size() << " approach levels." << std::endl;
for (size_type i=0; i < deltas.size(); i++) {
std::cout << "Solving approach level " << i+1 << " (" << deltas[i] << ")" << std::endl;
gmm::iteration iter(1e-6, 1, 40000);
getfem::default_newton_line_search ls;
block1.set_dirichlet_conditions_constant(md, -deltas[i]/2);
block2.set_dirichlet_conditions_constant(md, deltas[i]/2);
getfem::standard_solve(md, iter, getfem::rselect_linear_solver(md,"mumps"), ls);
std::stringstream ss_suffix;
ss_suffix << "_disp" << i;
block1.export_results_from_model(md, ss_suffix.str());
block2.export_results_from_model(md, ss_suffix.str());
}
return 0;
}
#include "getfem/getfem_regular_meshes.h"
#include "getfem/getfem_models.h"
#include "getfem/getfem_model_solvers.h"
#include <getfem/getfem_Coulomb_friction.h>
#include "getfem/getfem_export.h"
using bgeot::short_type;
using bgeot::size_type;
using bgeot::scalar_type;
using bgeot::base_small_vector;
typedef getfem::model_real_plain_vector plain_vector;
typedef getfem::model_real_sparse_matrix sparse_matrix;
typedef gmm::row_matrix<gmm::rsvector<scalar_type> > rsr_matrix;
void read_matrix_from_file(const std::string &fname, std::vector< std::vector<double> > &mat)
{
mat.clear();
std::ifstream inpfile(fname.c_str());
std::string strline;
while (std::getline(inpfile,strline)) {
unsigned int first = strline.find_first_not_of(" ");
if (first == std::string::npos ||
strline[first] == '#')
continue; // skip empty or commented out line
std::stringstream linestream(strline);
std::string strvalue;
std::vector<double> vecline(0);
while (std::getline(linestream, strvalue, ',')) {
double dvalue = std::atof(strvalue.c_str());
vecline.push_back(dvalue);
}
mat.push_back(vecline);
}
}
// mesh regions
enum { XP = 0,
XM = 1,
YP = 2,
YM = 3,
ZP = 4,
ZM = 5
};
class block
{
private:
int id;
double dx,dy,dz,lz;
double lambda,mu;
int nx,ny,nz,zsteps;
bool reversed;
double delta;
std::vector< std::vector<double> > roughnessmat;
std::vector<getfem::mesh> meshes;
std::vector<getfem::mesh_fem> mfus, mfds;
std::vector<getfem::mesh_im> mims;
public:
block(int id_, double dx_, double dy_, double dz_, double lz_, int nz_,
double E, double nu, const std::string &surface_file, bool uniform_dz=false)
: id(id_), dx(dx_), dy(dy_), dz(dz_), lz(lz_), nx(0), ny(0), nz(nz_),
reversed(false), delta(0),
meshes(0), mfus(0), mfds(0), mims(0)
{
zsteps = int(floor(log2(1+lz/(dz*nz))));
if (uniform_dz)
dz = lz/(nz*(pow(2,zsteps)-1));
lambda = E*nu / ((1+nu)*(1-2*nu));
mu = E/(2*(1+nu));
read_matrix_from_file(surface_file, roughnessmat);
meshes.resize(zsteps);
mfus.resize(zsteps);
mfds.resize(zsteps);
mims.resize(zsteps);
}
int nx_max(void) { return ny_max() >= 0 ? roughnessmat[0].size()-1 : 0; }
int ny_max(void) { return roughnessmat.size()-1; }
void build_meshes(int nx_, int ny_, int dir, int nnx=0, int nny=0)
{
nx = std::min(nx_, nx_max());
ny = std::min(ny_, ny_max());
if (nnx < 1 || nnx > nx)
nnx = 1;
if (nny < 1 || nny > ny)
nny = 1;
int nx0 = nx/nnx;
nx = nx0*nnx;
int ny0 = ny/nny;
ny = ny0*nny;
reversed = (dir < 0);
bgeot::base_node org;
std::vector<bgeot::base_small_vector> vect(3);
std::vector<int> ref(3);
for (int i=0; i < zsteps; i++) {
size_type nX = nnx * std::max(1,(int)round(nx0/pow(2.,i)));
size_type nY = nny * std::max(1,(int)round(ny0/pow(2.,i)));
ref[0] = nX; ref[1] = nY; ref[2] = nz;
double dX = (dx*nx)/nX;
double dY = (dy*ny)/nY;
double dZ = dz*pow(2,i);
double stretch_factor=1.;
if (i == zsteps-1) // for non-uniform dz
stretch_factor = (lz - (pow(2,i)-1)*dz*nz)/(dZ*nz);
vect[0] = bgeot::base_small_vector(dX, 0.0, 0.0);
vect[1] = bgeot::base_small_vector(0.0, dY, 0.0);
vect[2] = bgeot::base_small_vector(0.0, 0.0, stretch_factor*dZ);
if (dir >= 0)
org = bgeot::base_node(0.0, 0.0, (dZ-dz)*nz);
else
org = bgeot::base_node(0.0, 0.0, -((1+stretch_factor)*dZ-dz)*nz);
getfem::parallelepiped_regular_mesh(meshes[i], 3, org, vect.begin(), ref.begin());
// define regions on the boundary
getfem::mesh &m = meshes[i];
for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
int nb_faces = m.structure_of_convex(cv)->nb_faces();
for (short_type fc=0; fc < nb_faces; fc++) {
if (!m.is_convex_having_neighbour(cv,fc)) {
base_small_vector un = m.normal_of_face_of_convex(cv,fc,0);
gmm::scale(un, 1/gmm::vect_norm2(un));
if (gmm::abs(un[0] - 1.0) < 1.0E-3)
m.region(XP).add(cv,fc);
else if (gmm::abs(un[0] + 1.0) < 1.0E-3)
m.region(XM).add(cv,fc);
else if (gmm::abs(un[1] - 1.0) < 1.0E-3)
m.region(YP).add(cv,fc);
else if (gmm::abs(un[1] + 1.0) < 1.0E-3)
m.region(YM).add(cv,fc);
else if (gmm::abs(un[2] - 1.0) < 1.0E-3)
m.region(ZP).add(cv,fc);
else if (gmm::abs(un[2] + 1.0) < 1.0E-3)
m.region(ZM).add(cv,fc);
}
} // fc
} // cv
}
// apply surface roughness topography
for (dal::bv_visitor ip(meshes[0].points_index()); !ip.finished(); ++ip) {
for (int i=nz-1; i >= 0; i--) {
double z = ((dir >= 0) ? dz : -dz) * i;
double scale = ((dir >= 0) ? 1 : -1) * double(nz-i)/double(nz);
if (gmm::abs(meshes[0].points()[ip][2]-z) < 1e-10) {
size_type ix=round(meshes[0].points()[ip][0]/dx);
size_type iy=round(meshes[0].points()[ip][1]/dy);
meshes[0].points()[ip][2] -= scale * roughnessmat[iy][ix];
}
}
}
}
void define_mesh_fems(void)
{
// mesh_fem's and integration methods
for (int i=0; i < zsteps; i++) {
mfus[i].init_with_mesh(meshes[i]);
mfds[i].init_with_mesh(meshes[i]);
mfus[i].set_qdim(3);
mfds[i].set_qdim(1);
mfus[i].set_classical_finite_element(1);
mfds[i].set_classical_finite_element(1);
mims[i].init_with_mesh(meshes[i]);
mims[i].set_integration_method(getfem::int_method_descriptor("IM_GAUSS_PARALLELEPIPED(3,5)"));
}
}
void add_linear_elasticity_to_model(getfem::model &md)
{
md.add_initialized_scalar_data(lambda_str(), lambda);
md.add_initialized_scalar_data(mu_str(), mu);
for (int i=0; i < zsteps; i++) {
md.add_fem_variable(var_str(i), mfus[i]);
getfem::add_isotropic_linearized_elasticity_brick(md, mims[i], var_str(i), lambda_str(), mu_str());
}
}
void add_dirichlet_conditions_to_model(getfem::model &md, double disp)
{
delta = disp;
int i=zsteps-1;
int reg=(reversed ? ZM : ZP);
plain_vector W(mfus[i].nb_dof());
gmm::clear(W);
dal::bit_vector cn=mfus[i].dof_on_region(reg);
for (dal::bv_visitor dof(cn); !dof.finished(); ++dof)
if (dof % 3 == 2)
W[dof] = disp;
md.add_initialized_fem_data(dirichlet_displacements_str(), mfus[i], W);
getfem::add_Dirichlet_condition_with_multipliers
(md, mims[i], var_str(i), mfus[i], reg, dirichlet_displacements_str());
//getfem::add_Dirichlet_condition_with_penalization
// (md, mims[i], var_str(i), 1e9, reg, dirichlet_displacements_str());
}
void set_dirichlet_conditions_constant(getfem::model &md, double disp)
{
delta = disp;
int i=zsteps-1;
int reg=(reversed ? ZM : ZP);
plain_vector W(mfus[i].nb_dof());
gmm::clear(W);
dal::bit_vector cn=mfus[i].dof_on_region(reg);
for (dal::bv_visitor dof(cn); !dof.finished(); ++dof)
if (dof % 3 == 2)
W[dof] = disp;
gmm::copy(W, md.set_real_variable(dirichlet_displacements_str()));
}
void export_results_from_model(getfem::model &md, const std::string &suffix="")
{
for (int i=0; i < zsteps; i++) {
std::stringstream ss_id_i;
ss_id_i << id << "_" << i;
std::stringstream ss_fname;
ss_fname << "./output/roughcontact" << id << "_" << i << suffix << ".vtk";
getfem::vtk_export exp(ss_fname.str(), true);
exp.exporting(meshes[i]);
exp.write_point_data(mfus[i], md.real_variable(var_str(i)), "elastostatic_displacement" + ss_id_i.str());
plain_vector VonMisesStresses(mfds[i].nb_dof());
getfem::compute_isotropic_linearized_Von_Mises_or_Tresca
(md, var_str(i), lambda_str(), mu_str(), mfds[i], VonMisesStresses, false);
exp.write_point_data(mfds[i], VonMisesStresses, "von mises stresses");
}
}
const std::string lambda_str(void)
{
std::stringstream ss_lambda;
ss_lambda << "lambda" << id;
return ss_lambda.str();
}
const std::string mu_str(void)
{
std::stringstream ss_mu;
ss_mu << "mu" << id;
return ss_mu.str();
}
const std::string var_str(int i)
{
if (i >= 0 && i <= zsteps-1) {
std::stringstream ss_var;
ss_var << "u" << id << "_" << i;
return ss_var.str();
}
else
return "";
}
const std::string mortar_mult_str(int i)
{
if (i >= 1 && i <= zsteps-1) {
std::stringstream ss_mult;
ss_mult << "mortar_mult" << id << "_" << i;
return ss_mult.str();
}
else
return "";
}
const std::string periodicity_mult_str(void)
{
std::stringstream ss_mult;
ss_mult << "periodicity_mult" << id;
return ss_mult.str();
}
const std::string dirichlet_displacements_str()
{
std::stringstream ss_dirichlet(std::stringstream::out);
ss_dirichlet << "dirichlet_displacements" << id;
return ss_dirichlet.str();
}
const getfem::mesh_im& mim(int i)
{
return mims[i];
}
double suggested_augmentation_parameter(void)
{
return mu * (3*lambda + 2*mu) / (lambda + mu);
}
};
int main(int argc, char* argv[])
{
bgeot::md_param PARAM;
PARAM.read_command_line(argc, argv);
double dx = PARAM.real_value("DX","resolution in x direction [mm]");
double dy = PARAM.real_value("DY","resolution in y direction [mm]");
double lz1 = PARAM.real_value("LZ1","total height of the first body [mm]");
double lz2 = PARAM.real_value("LZ1","total height of the second body [mm]");
double dz1 = PARAM.real_value("DZMIN1","minimum element size in z direction (1) [mm]");
double dz2 = PARAM.real_value("DZMIN2","minimum element size in z direction (2) [mm]");
const std::vector<bgeot::md_param::param_value> &delta_params = PARAM.array_value("DELTA");
std::vector<double> deltas(delta_params.size());
for (size_type i=0; i < deltas.size(); ++i) {
GMM_ASSERT1(delta_params[i].type_of_param() == bgeot::md_param::REAL_VALUE,
"DELTA should be a scalar array.");
deltas[i] = delta_params[i].real();
}
int nz1 = PARAM.int_value("NZ_PER_LAYER1","number of elements per mesh layer (1)");
int nz2 = PARAM.int_value("NZ_PER_LAYER2","number of elements per mesh layer (2)");
double E1 = PARAM.real_value("E1","Young's Modulus (1) [N/mm²]");
double E2 = PARAM.real_value("E2","Young's Modulus (2) [N/mm²]");
double nu1 = PARAM.real_value("NU1","Poisson's ratio (1) [-]");
double nu2 = PARAM.real_value("NU2","Poisson's ratio (2) [-]");
std::string roughness1_file = PARAM.string_value("SURFACE1","filename for the 1st roughness topography");
std::string roughness2_file = PARAM.string_value("SURFACE2","filename for the 1st roughness topography");
block block1(1,dx,dy,dz1,lz1,nz1,E1,nu1,roughness1_file);
block block2(2,dx,dy,dz2,lz2,nz2,E2,nu2,roughness2_file);
size_type nx = std::min(block1.nx_max(), block2.nx_max());
size_type ny = std::min(block1.ny_max(), block2.ny_max());
int nnx = PARAM.int_value("PERIODIC_REPETITIONS_X");
int nny = PARAM.int_value("PERIODIC_REPETITIONS_Y");
// meshes
block1.build_meshes(nx,ny,1,nnx,nny);
block2.build_meshes(nx,ny,-1,nnx,nny);
// mesh_fem's
block1.define_mesh_fems();
block2.define_mesh_fems();
// model
getfem::model md;
block1.add_linear_elasticity_to_model(md);
block2.add_linear_elasticity_to_model(md);
block1.add_dirichlet_conditions_to_model(md,-deltas[0]/2);
block2.add_dirichlet_conditions_to_model(md,deltas[0]/2);
md.add_initialized_scalar_data("r", 1e6);
getfem::add_penalized_contact_between_nonmatching_meshes_brick
(md, block1.mim(0), block1.var_str(0), block2.var_str(0), "r", ZM, ZP);
std::cout << "Solving model with " << md.nb_dof() << " degrees of freedom at "
<< deltas.size() << " approach levels." << std::endl;
for (size_type i=0; i < deltas.size(); i++) {
std::cout << "Solving approach level " << i+1 << " (" << deltas[i] << ")" << std::endl;
gmm::iteration iter(1e-6, 1, 40000);
getfem::default_newton_line_search ls;
block1.set_dirichlet_conditions_constant(md, -deltas[i]/2);
block2.set_dirichlet_conditions_constant(md, deltas[i]/2);
getfem::standard_solve(md, iter, getfem::rselect_linear_solver(md,"mumps"), ls);
std::stringstream ss_suffix;
ss_suffix << "_disp" << i;
block1.export_results_from_model(md, ss_suffix.str());
block2.export_results_from_model(md, ss_suffix.str());
}
return 0;
}
_______________________________________________ Getfem-users mailing list [email protected] https://mail.gna.org/listinfo/getfem-users
