[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Yves Renard |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Fri, 15 Dec 2017 07:59:40 -0500 (EST) |
branch: master
commit ba95c87f531bf4187bdb2de05819715854515abf
Author: Yves Renard <address@hidden>
Date: Fri Dec 15 13:58:50 2017 +0100
Adding the possibility to print gfsparse matrices
---
interface/src/gf_spmat_get.cc | 12 +-
interface/tests/python/demo_dynamic_contact_1D.py | 327 ++++------------------
2 files changed, 66 insertions(+), 273 deletions(-)
diff --git a/interface/src/gf_spmat_get.cc b/interface/src/gf_spmat_get.cc
index 0815c32..fa5b664 100644
--- a/interface/src/gf_spmat_get.cc
+++ b/interface/src/gf_spmat_get.cc
@@ -357,9 +357,15 @@ void gf_spmat_get(getfemint::mexargs_in& m_in,
@*/
sub_command
("char", 0, 0, 0, 1,
- GMM_ASSERT1(false, "Sorry, function to be done");
- // std::string s = ...;
- // out.pop().from_string(s.c_str());
+ std::stringstream s;
+ if (gsp.storage() == getfemint::gsparse::WSCMAT) {
+ if (!gsp.is_complex()) s << gsp.wsc(scalar_type());
+ else s << gsp.wsc(complex_type());
+ } else {
+ if (!gsp.is_complex()) s << gsp.csc(scalar_type());
+ else s << gsp.csc(complex_type());
+ }
+ out.pop().from_string(s.str().c_str());
);
diff --git a/interface/tests/python/demo_dynamic_contact_1D.py
b/interface/tests/python/demo_dynamic_contact_1D.py
index 234e01c..ab06b0f 100644
--- a/interface/tests/python/demo_dynamic_contact_1D.py
+++ b/interface/tests/python/demo_dynamic_contact_1D.py
@@ -30,7 +30,6 @@ import getfem as gf
import numpy as np
-
# Parameters
NX=20 # Number of elements
T = 10 # Simulation time
@@ -47,292 +46,80 @@ theta = 1.0 # Theta-scheme parameter
singular_mass = 0 # singular mass or not
scheme = 0
+version = 0; # 0 = penalized contact ...
# Mesh
-m=gf.Mesh('cartesian', [0:1/NX:1]);
-
-
-exit(1);
-
-
-
-# Parameters of the model
-clambda = 1. # Lame coefficient
-cmu = 1. # Lame coefficient
-friction_coeff = 0.4 # coefficient of friction
-vertical_force = 0.05 # Volumic load in the vertical direction
-r = 10. # Augmentation parameter
-condition_type = 0 # 0 = Explicitely kill horizontal rigid displacements
- # 1 = Kill rigid displacements using a global penalization
- # 2 = Add a Dirichlet condition on the top of the structure
-penalty_parameter = 1E-6 # Penalization coefficient for the global
penalization
-
-if d == 2:
- cpoints = [0, 0] # constrained points for 2d
- cunitv = [1, 0] # corresponding constrained directions for 2d
-else:
- cpoints = [0, 0, 0, 0, 0, 0, 5, 0, 5] # constrained points for 3d
- cunitv = [1, 0, 0, 0, 1, 0, 0, 1, 0] # corresponding constrained
directions for 3d
-
-niter = 100 # Maximum number of iterations for Newton's algorithm.
-version = 13 # 1 : frictionless contact and the basic contact brick
- # 2 : contact with 'static' Coulomb friction and basic contact
brick
- # 3 : frictionless contact and the contact with a
- # rigid obstacle brick
- # 4 : contact with 'static' Coulomb friction and the contact
with a
- # rigid obstacle brick
- # 5 : frictionless contact and the integral brick
- # Newton and Alart-Curnier augmented lagrangian,
- # unsymmetric version
- # 6 : frictionless contact and the integral brick
- # Newton and Alart-Curnier augmented lagrangian, symmetric
- # version.
- # 7 : frictionless contact and the integral brick
- # Newton and Alart-Curnier augmented lagrangian,
- # unsymmetric version with an additional augmentation.
- # 8 : frictionless contact and the integral brick
- # New unsymmetric method.
- # 9 : frictionless contact and the integral brick : Uzawa
- # on the Lagrangian augmented by the penalization term.
- # 10 : contact with 'static' Coulomb friction and the integral
brick
- # Newton and Alart-Curnier augmented lagrangian,
- # unsymmetric version.
- # 11 : contact with 'static' Coulomb friction and the integral
brick
- # Newton and Alart-Curnier augmented lagrangian,
- # nearly symmetric version.
- # 12 : contact with 'static' Coulomb friction and the integral
brick
- # Newton and Alart-Curnier augmented lagrangian,
- # unsymmetric version with an additional augmentation.
- # 13 : contact with 'static' Coulomb friction and the integral
brick
- # New unsymmetric method.
- # 14 : contact with 'static' Coulomb friction and the integral
brick : Uzawa
- # on the Lagrangian augmented by the penalization term.
- # 15 : penalized contact with 'static' Coulomb friction (r is
the penalization
- # coefficient).
-
-# Signed distance representing the obstacle
-if d == 2:
- obstacle = 'y'
-else:
- obstacle = 'z'
+m=gf.Mesh('cartesian', np.arange(0,1+1./NX,1./NX))
# Selection of the contact and Dirichlet boundaries
-GAMMAC = 1
-GAMMAD = 2
+GAMMAC = 1; GAMMAD = 2
border = m.outer_faces()
normals = m.normal_of_faces(border)
-contact_boundary = border[:,np.nonzero(normals[d-1] < -0.01)[0]]
+contact_boundary = border[:,np.nonzero(normals[0] < -0.01)[0]]
m.set_region(GAMMAC, contact_boundary)
-contact_boundary = border[:,np.nonzero(normals[d-1] > 0.01)[0]]
+contact_boundary = border[:,np.nonzero(normals[0] > 0.01)[0]]
m.set_region(GAMMAD, contact_boundary)
# Finite element methods
-u_degree = 2
-lambda_degree = 2
-
-mfu = gf.MeshFem(m, d)
+mfu = gf.MeshFem(m)
mfu.set_classical_fem(u_degree)
mfd = gf.MeshFem(m, 1)
mfd.set_classical_fem(u_degree)
-mflambda = gf.MeshFem(m, 1) # used only by version 5 to 13
-mflambda.set_classical_fem(lambda_degree)
-
-mfvm = gf.MeshFem(m, 1)
-mfvm.set_classical_discontinuous_fem(u_degree-1)
-
# Integration method
mim = gf.MeshIm(m, 4)
-if d == 2:
- mim_friction = gf.MeshIm(m,
- gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(4),4)'))
-else:
- mim_friction = gf.MeshIm(m,
- gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TETRAHEDRON(5),4)'))
-
-# Volumic density of force
-nbdofd = mfd.nbdof()
-nbdofu = mfu.nbdof()
-F = np.zeros(nbdofd*d)
-F[d-1:nbdofd*d:d] = -vertical_force;
-
-# Elasticity model
-md = gf.Model('real')
-md.add_fem_variable('u', mfu)
-md.add_initialized_data('cmu', [cmu])
-md.add_initialized_data('clambda', [clambda])
-md.add_isotropic_linearized_elasticity_brick(mim, 'u', 'clambda', 'cmu')
-md.add_initialized_fem_data('volumicload', mfd, F)
-md.add_source_term_brick(mim, 'u', 'volumicload')
-
-if condition_type == 2:
- Ddata = np.zeros(d)
- Ddata[d-1] = -5
- md.add_initialized_data('Ddata', Ddata)
- md.add_Dirichlet_condition_with_multipliers(mim, 'u', u_degree, GAMMAD,
'Ddata')
-elif condition_type == 0:
- md.add_initialized_data('cpoints', cpoints)
- md.add_initialized_data('cunitv', cunitv)
- md.add_pointwise_constraints_with_multipliers('u', 'cpoints', 'cunitv')
-elif condition_type == 1:
- # Small penalty term to avoid rigid motion (should be replaced by an
- # explicit treatment of the rigid motion with a constraint matrix)
- md.add_initialized_data('penalty_param', [penalty_parameter])
- md.add_mass_brick(mim, 'u', 'penalty_param')
-
-# The contact condition
-
-cdof = mfu.dof_on_region(GAMMAC)
-nbc = cdof.shape[0] / d
-
-solved = False
-if version == 1 or version == 2: # defining the matrices BN and BT by hand
- contact_dof = cdof[d-1:nbc*d:d]
- contact_nodes = mfu.basic_dof_nodes(contact_dof)
- BN = gf.Spmat('empty', nbc, nbdofu)
- ngap = np.zeros(nbc)
- for i in range(nbc):
- BN[i, contact_dof[i]] = -1.
- ngap[i] = contact_nodes[d-1, i]
-
- if version == 2:
- BT = gf.Spmat('empty', nbc*(d-1), nbdofu)
- for i in range(nbc):
- for j in range(d-1):
- BT[j+i*(d-1), contact_dof[i]-d+j+1] = 1.0
-
- md.add_variable('lambda_n', nbc)
- md.add_initialized_data('r', [r])
- if version == 2:
- md.add_variable('lambda_t', nbc * (d-1))
- md.add_initialized_data('friction_coeff', [friction_coeff])
-
- md.add_initialized_data('ngap', ngap)
- md.add_initialized_data('alpha', np.ones(nbc))
- if version == 1:
- md.add_basic_contact_brick('u', 'lambda_n', 'r', BN, 'ngap', 'alpha', 1)
- else:
- md.add_basic_contact_brick('u', 'lambda_n', 'lambda_t', 'r', BN, BT,
'friction_coeff', 'ngap', 'alpha', 1);
-
-elif version == 3 or version == 4: # BN and BT defined by the contact brick
-
- md.add_variable('lambda_n', nbc)
- md.add_initialized_data('r', [r])
- if version == 3:
- md.add_nodal_contact_with_rigid_obstacle_brick(mim, 'u', 'lambda_n',
'r', GAMMAC, obstacle, 1);
- else:
- md.add_variable('lambda_t', nbc*(d-1))
- md.add_initialized_data('friction_coeff', [friction_coeff])
- md.add_nodal_contact_with_rigid_obstacle_brick(mim, 'u', 'lambda_n',
'lambda_t', 'r',
- 'friction_coeff', GAMMAC,
obstacle, 1)
-
-elif version >= 5 and version <= 8: # The integral version, Newton
-
- ldof = mflambda.dof_on_region(GAMMAC)
- mflambda_partial = gf.MeshFem('partial', mflambda, ldof)
- md.add_fem_variable('lambda_n', mflambda_partial)
- md.add_initialized_data('r', [r])
- OBS = mfd.eval(obstacle)
- md.add_initialized_fem_data('obstacle', mfd, OBS)
- md.add_integral_contact_with_rigid_obstacle_brick(mim_friction, 'u',
'lambda_n',
- 'obstacle', 'r', GAMMAC,
version-4);
-
-elif version == 9: # The integral version, Uzawa on the augmented Lagrangian
-
- ldof = mflambda.dof_on_region(GAMMAC)
- mflambda_partial = gf.MeshFem('partial', mflambda, ldof)
- nbc = mflambda_partial.nbdof()
- M = gf.asm_mass_matrix(mim, mflambda_partial, mflambda_partial, GAMMAC)
- lambda_n = np.zeros(nbc)
- md.add_initialized_fem_data('lambda_n', mflambda_partial, lambda_n)
- md.add_initialized_data('r', [r])
- OBS = mfd.eval(obstacle) # np.array([mfd.eval(obstacle)])
- md.add_initialized_fem_data('obstacle', mfd, OBS)
- md.add_penalized_contact_with_rigid_obstacle_brick \
- (mim_friction, 'u', 'obstacle', 'r', GAMMAC, 2, 'lambda_n')
-
- for ii in range(100):
- print('iteration %d' % (ii+1))
- md.solve('max_res', 1E-9, 'max_iter', niter)
- U = md.get('variable', 'u')
- lambda_n_old = lambda_n
- sol = gf.linsolve_superlu(M,
gf.asm_integral_contact_Uzawa_projection(GAMMAC, mim_friction, mfu, U,
mflambda_partial, lambda_n, mfd, OBS, r))
- lambda_n = sol[0].transpose()
- md.set_variable('lambda_n', lambda_n)
- difff = max(abs(lambda_n-lambda_n_old))[0]/max(abs(lambda_n))[0]
- print('diff : %g' % difff)
- if difff < penalty_parameter:
- break
-
- solved = True
-
-elif version >= 10 and version <= 13: # The integral version with friction,
Newton
-
- mflambda.set_qdim(d);
- ldof = mflambda.dof_on_region(GAMMAC)
- mflambda_partial = gf.MeshFem('partial', mflambda, ldof)
- md.add_fem_variable('lambda', mflambda_partial)
- md.add_initialized_data('r', [r])
- md.add_initialized_data('friction_coeff', [friction_coeff])
- OBS = mfd.eval(obstacle)
- md.add_initialized_fem_data('obstacle', mfd, OBS)
- md.add_integral_contact_with_rigid_obstacle_brick \
- (mim_friction, 'u', 'lambda', 'obstacle', 'r', 'friction_coeff', GAMMAC,
version-9)
-
-elif version == 14: # The integral version, Uzawa on the augmented Lagrangian
with friction
-
- mflambda.set_qdim(d)
- ldof = mflambda.dof_on_region(GAMMAC)
- mflambda_partial = gf.MeshFem('partial', mflambda, ldof)
- nbc = mflambda_partial.nbdof()
- md.add_initialized_data('friction_coeff', [friction_coeff])
- M = gf.asm_mass_matrix(mim, mflambda_partial, mflambda_partial, GAMMAC)
- lambda_nt = np.zeros(nbc)
- md.add_initialized_fem_data('lambda', mflambda_partial, lambda_nt)
- md.add_initialized_data('r', [r])
- OBS = mfd.eval(obstacle)
- md.add_initialized_fem_data('obstacle', mfd, OBS)
- md.add_penalized_contact_with_rigid_obstacle_brick \
- (mim_friction, 'u', 'obstacle', 'r', 'friction_coeff', GAMMAC, 2,
'lambda')
-
- for ii in range(100):
- print('iteration %d' % (ii+1))
- md.solve('max_res', 1E-9, 'max_iter', niter)
- U = md.get('variable', 'u')
- lambda_nt_old = lambda_nt
- sol = gf.linsolve_superlu(M,
- gf.asm_integral_contact_Uzawa_projection(
- GAMMAC, mim_friction, mfu, U, mflambda_partial, lambda_nt, mfd, OBS,
r, friction_coeff))
- lambda_nt = sol[0].transpose()
- md.set_variable('lambda', lambda_nt)
- difff = max(abs(lambda_nt-lambda_nt_old))[0]/max(abs(lambda_nt))[0]
- print('diff : %g' % difff)
- if difff < penalty_parameter:
- break
-
- solved = True
-
-elif version == 15:
-
- md.add_initialized_data('r', [r])
- md.add_initialized_data('friction_coeff', [friction_coeff])
- OBS = mfd.eval(obstacle)
- md.add_initialized_fem_data('obstacle', mfd, OBS);
- md.add_penalized_contact_with_rigid_obstacle_brick \
- (mim_friction, 'u', 'obstacle', 'r', 'friction_coeff', GAMMAC)
-
-else:
- print('Inexistent version')
-
-# Solve the problem
-if not solved:
- md.solve('max_res', 1E-9, 'very noisy', 'max_iter', niter, 'lsearch',
'default') #, 'with pseudo potential')
-
-U = md.get('variable', 'u')
-# LAMBDA = md.get('variable', 'lambda_n')
-VM = md.compute_isotropic_linearized_Von_Mises_or_Tresca('u', 'clambda',
'cmu', mfvm)
-mfd.export_to_vtk('static_contact.vtk', 'ascii', mfvm, VM, 'Von Mises
Stress', mfu, U, 'Displacement')
+# Mass and stiffness matrices
+md = gf.Model('real'); md.add_fem_variable('u', mfu);
+M = gf.asm_generic(mim, 2, 'u*Test_u', -1, md)
+K = gf.asm_generic(mim, 2, 'Grad_u*Grad_Test_u', -1, md)
+
+# Dirichlet condition on the top
+for i in range(0, u_degree+1): K[0,i] = 0.;
+for i in range(0, u_degree+1): M[0,i] = 0.;
+M[0,0] = K[0,0] = 1.;
+
+print uExact(0,0)
+
+
+print K
+print K[1,1]
+
+
+
+
+def uExact(x, t):
+ # The solution is periodic of period 3
+ # Shift the time 't' with t=0 : beginning of the period
+ tp = rem(t,3.)
+ # The solution has 3 phases
+ # Shift the time 'tp' with t=0 : beginning of a phase
+ # and get also the phase number
+ tf = rem(tp,1.)
+ nf = floor(tp)
+ # Get the index of the zone in each phase : I, II, III, IV
+ # (zones are given according to characteristics of the wave equation)
+ if (tf<=x):
+ if (tf<=(1-x)): zone = 1; else: zone = 3;
+ else:
+ if (tf<=(1-x)): zone = 2; else: zone = 4;
+ # Solution according to the Phase (1,2,3) and the zone (I, II, III, IV)
+ if nf == 0:
+ if zone == 1: u = 1./2.-x/2.; dxu = -1./2.; dtu = 0.;
+ elif zone == 2: u = 1./2.-tf/2.; dxu = 0.; dtu = -1./2.;
+ elif zone == 3: u = 1./2.-x/2.; dxu = -1./2.; dtu = 0;
+ elif zone == 4: u = 1./2.-tf/2.; dxu = 0; dtu = -1./2.;
+ elif nf == 1:
+ if zone == 1: u = -tf/2.; dxu = 0; dtu = -1./2.
+ elif zone == 2: u = -x/2.; dxu = -1./2.; dtu = 0;
+ elif zone == 3: u = -1./2.+x/2.; dxu = 1./2.; dtu = 0;
+ elif zone == 4: u = -1./2.+tf/2.; dxu = 0; dtu = 1./2.;
+ elif nf == 2:
+ if zone == 1: u = tf/2.; dxu = 0; dtu = 1./2.;
+ elif zone == 2: u = tf/2.; dxu = 0; dtu = 1./2.;
+ elif zone == 3: u = 1./2.-x/2.; dxu = -1./2.; dtu = 0;
+ elif zone == 4: u = 1./2.-x/2.; dxu = -1./2.; dtu = 0.
+ return u