You are not logged in. Please login or register.
OOFEM forum → Developers Discussion → Contact element implementation
You should be looking at the base classes, SparseMtrx:
answer->assemble(rloc, cloc, tangent);
is how you write to it.
Your code shouldn't try to handle anything else. The matrix is zeroed by the solver between each iteration.
You should just assemble your contributions to the global stiffness, essentially like how you have right now in the code.
You should be looking at the base classes, SparseMtrx:
answer->assemble(rloc, cloc, tangent);
is how you write to it.Your code shouldn't try to handle anything else. The matrix is zeroed by the solver between each iteration.
You should just assemble your contributions to the global stiffness, essentially like how you have right now in the code.
That's what I'm doing, and it can be seen from the code, but the implementation is taken for some reason from skyline.c. Also if I print the matrix, it gives the full structure, however I cannot access to the some elements I have OOFEM_ERROR("you're asking the element which is not presented in the matrix (something like this))... I'll have a look into it tomorrow.
The skyline matrix type doesnt allow for dynamic resizing, sobit must be preallocated. Im typing this on a tablet right now so i cant easily check, but you should be overloading the method for requesting the location arrays in you b.c.
The skyline matrix type doesnt allow for dynamic resizing, sobit must be preallocated. Im typing this on a tablet right now so i cant easily check, but you should be overloading the method for requesting the location arrays in you b.c.
I checked it, and it seems that problem is because in the constructor of the staticstructural the sparseMtrxType is set to sparseMtrxType(SMT_Skyline).
StaticStructural :: StaticStructural(int i, EngngModel *_master) : StructuralEngngModel(i, _master),
internalForces(),
eNorm(),
sparseMtrxType(SMT_Skyline)
{
ndomains = 1;
solverType = 0;
mRecomputeStepAfterPropagation = false;
}I found the types of the SparseMtrxType. What is the correct type? I believe I don't allow to change anything in the staticstructural.c ...
/**
* Enumerative type used to identify the sparse matrix type.
*/
enum SparseMtrxType {
SMT_Skyline, ///< Symmetric skyline.
SMT_SkylineU, ///< Unsymmetric skyline.
SMT_CompCol, ///< Compressed column.
SMT_DynCompCol, ///< Dynamically growing compressed column.
SMT_SymCompCol, ///< Symmetric compressed column.
SMT_DynCompRow, ///< Dynamically growing compressed row.
SMT_SpoolesMtrx, ///< Spooles sparse mtrx representation.
SMT_PetscMtrx, ///< PETSc library mtrx representation.
SMT_DSS_sym_LDL, ///< Richard Vondracek's sparse direct solver.
SMT_DSS_sym_LL, ///< Richard Vondracek's sparse direct solver.
SMT_DSS_unsym_LU ///< Richard Vondracek's sparse direct solver.
};What is "sobit"? I request the location array in the assemble function itself, so the call answer.assemble(r_loc, c_loc, Ke); causes the problem since it tries to add elements which is not specified.
Regards,
Vitalii
Phone typo, "sobit" ---> "so it"
No.. there is no problem with skyline.
They should all work. There are methods you overload for your boundary conditions that the matrices use to determine what they should be preallocating.
It's vitally important that you never rename, or change *anything* in the function prototypes of the methods you overload. They must be exactly as specified
virtual void giveLocationArrays(std :: vector< IntArray > &rows, std :: vector< IntArray > &cols, CharType type,
const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s) { }I pointed this out last a month ago:
For this to work with OOFEM, it's not enough to just name the methods correctly. They must be exactly as defined
void NTS :: giveLocationArrays(std :: vector< IntArray > &rows, std :: vector< IntArray > &cols, CharType type, const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s)Whether or not this is the cause, I don't know. Perhaps the matrix type you use can't handle non-preallocated storage (which is what the "giveLocationArrays" method is used for)
and it still very much applies.
Hi Mikael!
It finally does work! Thanks a lot. Convergence is OK with the tangent implementation.
Total number of solution steps 11
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 1, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: *1.668e-013 D_v: *1.636e+007
NRSolver: 1 D_u: 4.591e-003 D_v: 5.018e-003
NRSolver: 2 D_u: 4.143e-008 D_v: 2.017e-008
EngngModel info: user time consumed by solution step 1: 0.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 2, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 4.738e-001 D_v: 2.093e-001
NRSolver: 1 D_u: 2.647e-001 D_v: 7.053e-002
NRSolver: 2 D_u: 3.362e-001 D_v: 3.200e-001
NRSolver: 3 D_u: 2.391e-003 D_v: 1.846e-003
NRSolver: 4 D_u: 5.023e-008 D_v: 2.659e-008
EngngModel info: user time consumed by solution step 2: 0.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 3, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 2.260e-002 D_v: 1.463e-002
NRSolver: 1 D_u: 7.044e-003 D_v: 4.626e-003
NRSolver: 2 D_u: 8.257e-005 D_v: 6.765e-006
EngngModel info: user time consumed by solution step 3: 0.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 4, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 2.424e-002 D_v: 4.001e-002
NRSolver: 1 D_u: 1.450e-002 D_v: 1.672e-002
NRSolver: 2 D_u: 1.010e-003 D_v: 2.170e-004
NRSolver: 3 D_u: 7.789e-007 D_v: 1.531e-007
EngngModel info: user time consumed by solution step 4: 0.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 5, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 5.917e-003 D_v: 4.070e-003
NRSolver: 1 D_u: 1.640e-004 D_v: 5.068e-006
EngngModel info: user time consumed by solution step 5: 0.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 6, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 4.712e-003 D_v: 3.497e-003
NRSolver: 1 D_u: 9.571e-005 D_v: 1.754e-005
EngngModel info: user time consumed by solution step 6: 0.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 7, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 4.069e-003 D_v: 3.099e-003
NRSolver: 1 D_u: 6.833e-005 D_v: 2.609e-005
EngngModel info: user time consumed by solution step 7: 0.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 8, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 3.575e-003 D_v: 2.782e-003
NRSolver: 1 D_u: 5.414e-005 D_v: 2.909e-005
EngngModel info: user time consumed by solution step 8: 0.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 9, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 3.192e-003 D_v: 2.526e-003
NRSolver: 1 D_u: 4.465e-005 D_v: 3.031e-005
EngngModel info: user time consumed by solution step 9: 0.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 10, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 2.891e-003 D_v: 2.315e-003
NRSolver: 1 D_u: 3.741e-005 D_v: 3.093e-005
EngngModel info: user time consumed by solution step 10: 0.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 11, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 2.649e-003 D_v: 2.139e-003
NRSolver: 1 D_u: 3.158e-005 D_v: 3.133e-005
EngngModel info: user time consumed by solution step 11: 0.00sThe video with the results are in the attachment.
My next step is the implementation of the contact element with friction. There I'll have the non-symmetrical tangent matrix, will it pose the problem? skyline.c stores only upper triangle and staticstructural has it in the constructor. Is it possible to work with skylineu.c for instance?
P.S. It's not straightforward that giveLocationArray uses to make the internal structure of the sparse matrix, now I finally understood it.
Regards and many thanks again,
Vitalii
Huh? "giveLocationArray" is *used by* the sparse matrices to construct the internal structure, not the other way around.
Also, the documentation does say
/**
* Gives a list of location arrays that will be assembled.
* This should only be used to construct zero structure in sparse matrices.
* The rows and columns location arrays returned in tuples (stored in vector),
* allowing to efficiently assemble and allocate off-diagonal blocks.
* The nonzero entries are assembled and allocated for entries at (rows[i], cols[i]) positions.
* @param rows List of location arrays for r_s.
* @param cols List of location arrays for c_s.
* @param type Type of matrix to assemble.
* @param r_s Row numbering scheme.
* @param c_s Column numbering scheme.
*/Skyline is only the default matrix type.
http://www.oofem.org/resources/doc/oofe … ode26.html
Huh? "giveLocationArray" is *used by* the sparse matrices to construct the internal structure, not the other way around.
Also, the documentation does say/** * Gives a list of location arrays that will be assembled. * This should only be used to construct zero structure in sparse matrices. * The rows and columns location arrays returned in tuples (stored in vector), * allowing to efficiently assemble and allocate off-diagonal blocks. * The nonzero entries are assembled and allocated for entries at (rows[i], cols[i]) positions. * @param rows List of location arrays for r_s. * @param cols List of location arrays for c_s. * @param type Type of matrix to assemble. * @param r_s Row numbering scheme. * @param c_s Column numbering scheme. */Skyline is only the default matrix type.
http://www.oofem.org/resources/doc/oofe … ode26.html
OK 
, sadly there is no clear road map for the contact element implementation. 
I'll come back with the friction element implementation, hopefully without any further questions.
Will it be interesting to add the 2D-contact elements into OOFEM? Of course, I'll do some tests to verify the elements. There are not so many open-source code with this kind of functionallity...
There is really just three parts.
1. Assemble the residual
2. Assemble the tangent
3. + give the location arrays for a dry assembly/preallocation of the sparse matrix.
which is very standard stuff.
I pointed you to SurfaceTensionBoundaryCondition, because, well.. it does all of this, and pretty much only this, but with different stiffness contributions.
I think there are a few things that need to be sorted out before it's ready for inclusion:
1. It should use sets correctly (using element boundaries, not lists of nodes where the order is assumed).
nodetoedgecontact 3 loadtimefunction 1 masterset 1 slaveset 2 k 3500
...
# Just making up some sets right now, but it should be element boundaries:
Set 1 elementboundaries 4 5 1 2 3 # element 5, edge 1, element 2 edge 3
Set 2 elementboundaries 4 7 2 4 1 # element 7, edge 2, element 4 edge 12. It should use the element for the basis function on the boundary. Again, this is what surfacetensionbc.C does:
FEInterpolation *fei = e->giveInterpolation();
...
fei->boundaryEvalN( N, side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );3. Also, coding style needs to be fixed.
4. It should use the octree spatial localizer for lookup. This is part is not critical to start with, but should be done eventually. (might need change the localizer to support large deformations here, but it's a solveable issue)
The thing is, I think more about what it will be in real usecases. You typically don't want to define the contact positions beforehand, you want it automatically looked up. So it's not a contact element. It's many-to-many contact condition, with many nodes and edges.
Hello Mikael,
I've finished the implementation of the Node-To-Segment element with friction. However, I have a problem with the convergence. If the problem doesn't have friction it converges as usual. I've spent quite some time on checking of the mathematics it should be correct. With friction I have the following print out of the solution:
StaticStructural :: solveYourselfAt - Solving step 14, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
FloatArray of size : 8
-1.768e-002 -5.621e-002 1.768e-002 -5.621e-002 2.621e-002 -4.169e-002 -2.621e-002 -4.169e-002
NRSolver: 0 D_u: 3.200e-001 D_v: 4.733e-004
NRSolver: 1 D_u: 1.272e-001 D_v: 1.944e-005
NRSolver: 2 D_u: 5.425e-001 D_v: 3.281e-005
NRSolver: 3 D_u: 1.625e-002 D_v: 5.230e-004
NRSolver: 4 D_u: 4.699e-001 D_v: 1.981e-007
NRSolver: 5 D_u: 3.125e-002 D_v: 1.291e-006
NRSolver: 6 D_u: 6.125e-001 D_v: 1.302e-006
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 7 D_u: 3.023e-001 D_v: 1.250e-004
NRSolver: 8 D_u: 3.960e-003 D_v: 1.255e-004
NRSolver: 9 D_u: 2.242e-001 D_v: 1.791e-008
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 10 D_u: 1.585e-001 D_v: 1.008e-007
NRSolver: 11 D_u: 1.183e-001 D_v: 1.878e-007
NRSolver: 12 D_u: 1.673e-001 D_v: 1.056e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 13 D_u: 2.086e-001 D_v: 1.061e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 14 D_u: 1.924e-001 D_v: 1.365e-007
NRSolver: 15 D_u: 1.100e-001 D_v: 2.856e-007
NRSolver: 16 D_u: 2.028e-001 D_v: 9.403e-008
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 17 D_u: 2.173e-001 D_v: 1.282e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 18 D_u: 2.154e-001 D_v: 1.578e-007
NRSolver: 19 D_u: 1.144e-001 D_v: 3.709e-007
NRSolver: 20 D_u: 2.266e-001 D_v: 1.055e-007
EngngModel info: user time consumed by solution step 14: 0.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 15, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
FloatArray of size : 8
-1.910e-002 -5.954e-002 1.910e-002 -5.954e-002 2.806e-002 -4.471e-002 -2.806e-002 -4.471e-002
NRSolver: 0 D_u: 3.210e-001 D_v: 4.424e-004
NRSolver: 1 D_u: 1.264e-001 D_v: 1.971e-005
NRSolver: 2 D_u: 5.405e-001 D_v: 3.471e-005
NRSolver: 3 D_u: 1.723e-002 D_v: 5.586e-004
NRSolver: 4 D_u: 4.773e-001 D_v: 2.350e-007
NRSolver: 5 D_u: 3.179e-002 D_v: 1.438e-006
NRSolver: 6 D_u: 6.102e-001 D_v: 1.452e-006
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 7 D_u: 3.014e-001 D_v: 1.333e-004
NRSolver: 8 D_u: 4.195e-003 D_v: 1.339e-004
NRSolver: 9 D_u: 2.206e-001 D_v: 2.133e-008
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 10 D_u: 1.593e-001 D_v: 1.061e-007
NRSolver: 11 D_u: 1.164e-001 D_v: 2.043e-007
NRSolver: 12 D_u: 1.681e-001 D_v: 1.100e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 13 D_u: 2.071e-001 D_v: 1.130e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 14 D_u: 1.921e-001 D_v: 1.450e-007
NRSolver: 15 D_u: 1.092e-001 D_v: 3.065e-007
NRSolver: 16 D_u: 2.025e-001 D_v: 9.968e-008
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 17 D_u: 2.163e-001 D_v: 1.369e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 18 D_u: 2.148e-001 D_v: 1.685e-007
NRSolver: 19 D_u: 1.139e-001 D_v: 3.968e-007
NRSolver: 20 D_u: 2.260e-001 D_v: 1.125e-007
EngngModel info: user time consumed by solution step 15: 0.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 16, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
FloatArray of size : 8
-2.052e-002 -6.283e-002 2.052e-002 -6.283e-002 2.988e-002 -4.771e-002 -2.988e-002 -4.771e-002
NRSolver: 0 D_u: 3.218e-001 D_v: 4.158e-004
NRSolver: 1 D_u: 1.257e-001 D_v: 1.994e-005
NRSolver: 2 D_u: 5.385e-001 D_v: 3.660e-005
NRSolver: 3 D_u: 1.818e-002 D_v: 5.936e-004
NRSolver: 4 D_u: 4.842e-001 D_v: 2.744e-007
NRSolver: 5 D_u: 3.230e-002 D_v: 1.589e-006
NRSolver: 6 D_u: 6.078e-001 D_v: 1.606e-006
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 7 D_u: 3.004e-001 D_v: 1.415e-004
NRSolver: 8 D_u: 4.422e-003 D_v: 1.422e-004
NRSolver: 9 D_u: 2.171e-001 D_v: 2.500e-008
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 10 D_u: 1.601e-001 D_v: 1.113e-007
NRSolver: 11 D_u: 1.146e-001 D_v: 2.209e-007
NRSolver: 12 D_u: 1.690e-001 D_v: 1.141e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 13 D_u: 2.055e-001 D_v: 1.198e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 14 D_u: 1.918e-001 D_v: 1.533e-007
NRSolver: 15 D_u: 1.084e-001 D_v: 3.272e-007
NRSolver: 16 D_u: 2.022e-001 D_v: 1.052e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 17 D_u: 2.154e-001 D_v: 1.456e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 18 D_u: 2.141e-001 D_v: 1.789e-007
NRSolver: 19 D_u: 1.134e-001 D_v: 4.225e-007
NRSolver: 20 D_u: 2.254e-001 D_v: 1.193e-007
EngngModel info: user time consumed by solution step 16: 0.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 17, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
FloatArray of size : 8
-2.195e-002 -6.609e-002 2.195e-002 -6.609e-002 3.168e-002 -5.067e-002 -3.168e-002 -5.067e-002
NRSolver: 0 D_u: 3.222e-001 D_v: 3.926e-004
NRSolver: 1 D_u: 1.251e-001 D_v: 2.015e-005
NRSolver: 2 D_u: 5.364e-001 D_v: 3.846e-005
NRSolver: 3 D_u: 1.909e-002 D_v: 6.282e-004
NRSolver: 4 D_u: 4.906e-001 D_v: 3.161e-007
NRSolver: 5 D_u: 3.278e-002 D_v: 1.745e-006
NRSolver: 6 D_u: 6.054e-001 D_v: 1.765e-006
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 7 D_u: 2.994e-001 D_v: 1.496e-004
NRSolver: 8 D_u: 4.639e-003 D_v: 1.503e-004
NRSolver: 9 D_u: 2.137e-001 D_v: 2.892e-008
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 10 D_u: 1.607e-001 D_v: 1.164e-007
NRSolver: 11 D_u: 1.128e-001 D_v: 2.377e-007
NRSolver: 12 D_u: 1.697e-001 D_v: 1.180e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 13 D_u: 2.040e-001 D_v: 1.265e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 14 D_u: 1.915e-001 D_v: 1.615e-007
NRSolver: 15 D_u: 1.076e-001 D_v: 3.478e-007
NRSolver: 16 D_u: 2.019e-001 D_v: 1.105e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 17 D_u: 2.144e-001 D_v: 1.542e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 18 D_u: 2.135e-001 D_v: 1.893e-007
NRSolver: 19 D_u: 1.129e-001 D_v: 4.478e-007
NRSolver: 20 D_u: 2.247e-001 D_v: 1.260e-007
EngngModel info: user time consumed by solution step 17: 0.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 18, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
FloatArray of size : 8
-2.338e-002 -6.931e-002 2.338e-002 -6.931e-002 3.346e-002 -5.361e-002 -3.346e-002 -5.361e-002
NRSolver: 0 D_u: 3.225e-001 D_v: 3.722e-004
NRSolver: 1 D_u: 1.245e-001 D_v: 2.033e-005
NRSolver: 2 D_u: 5.343e-001 D_v: 4.032e-005
NRSolver: 3 D_u: 1.996e-002 D_v: 6.623e-004
NRSolver: 4 D_u: 4.965e-001 D_v: 3.599e-007
NRSolver: 5 D_u: 3.322e-002 D_v: 1.905e-006
NRSolver: 6 D_u: 6.030e-001 D_v: 1.928e-006
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 7 D_u: 2.984e-001 D_v: 1.576e-004
NRSolver: 8 D_u: 4.847e-003 D_v: 1.583e-004
NRSolver: 9 D_u: 2.104e-001 D_v: 3.305e-008
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 10 D_u: 1.613e-001 D_v: 1.214e-007
NRSolver: 11 D_u: 1.111e-001 D_v: 2.547e-007
NRSolver: 12 D_u: 1.703e-001 D_v: 1.216e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 13 D_u: 2.026e-001 D_v: 1.332e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 14 D_u: 1.911e-001 D_v: 1.695e-007
NRSolver: 15 D_u: 1.068e-001 D_v: 3.682e-007
NRSolver: 16 D_u: 2.015e-001 D_v: 1.157e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 17 D_u: 2.134e-001 D_v: 1.626e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 18 D_u: 2.128e-001 D_v: 1.994e-007
NRSolver: 19 D_u: 1.124e-001 D_v: 4.729e-007
NRSolver: 20 D_u: 2.240e-001 D_v: 1.327e-007
EngngModel info: user time consumed by solution step 18: 1.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 19, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
FloatArray of size : 8
-2.481e-002 -7.248e-002 2.481e-002 -7.248e-002 3.521e-002 -5.650e-002 -3.521e-002 -5.650e-002
NRSolver: 0 D_u: 3.226e-001 D_v: 3.542e-004
NRSolver: 1 D_u: 1.240e-001 D_v: 2.049e-005
NRSolver: 2 D_u: 5.322e-001 D_v: 4.216e-005
NRSolver: 3 D_u: 2.079e-002 D_v: 6.960e-004
NRSolver: 4 D_u: 5.019e-001 D_v: 4.054e-007
NRSolver: 5 D_u: 3.363e-002 D_v: 2.068e-006
NRSolver: 6 D_u: 6.005e-001 D_v: 2.095e-006
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 7 D_u: 2.973e-001 D_v: 1.655e-004
NRSolver: 8 D_u: 5.046e-003 D_v: 1.662e-004
NRSolver: 9 D_u: 2.072e-001 D_v: 3.737e-008
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 10 D_u: 1.619e-001 D_v: 1.262e-007
NRSolver: 11 D_u: 1.094e-001 D_v: 2.718e-007
NRSolver: 12 D_u: 1.709e-001 D_v: 1.249e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 13 D_u: 2.011e-001 D_v: 1.397e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 14 D_u: 1.907e-001 D_v: 1.774e-007
NRSolver: 15 D_u: 1.061e-001 D_v: 3.885e-007
NRSolver: 16 D_u: 2.011e-001 D_v: 1.208e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 17 D_u: 2.124e-001 D_v: 1.710e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 18 D_u: 2.121e-001 D_v: 2.095e-007
NRSolver: 19 D_u: 1.119e-001 D_v: 4.977e-007
NRSolver: 20 D_u: 2.233e-001 D_v: 1.392e-007
EngngModel info: user time consumed by solution step 19: 0.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 20, metastep 1, (neq = 8)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
FloatArray of size : 8
-2.626e-002 -7.562e-002 2.626e-002 -7.562e-002 3.693e-002 -5.937e-002 -3.693e-002 -5.937e-002
NRSolver: 0 D_u: 3.225e-001 D_v: 3.383e-004
NRSolver: 1 D_u: 1.235e-001 D_v: 2.064e-005
NRSolver: 2 D_u: 5.301e-001 D_v: 4.399e-005
NRSolver: 3 D_u: 2.159e-002 D_v: 7.292e-004
NRSolver: 4 D_u: 5.068e-001 D_v: 4.525e-007
NRSolver: 5 D_u: 3.400e-002 D_v: 2.235e-006
NRSolver: 6 D_u: 5.979e-001 D_v: 2.266e-006
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 7 D_u: 2.962e-001 D_v: 1.732e-004
NRSolver: 8 D_u: 5.236e-003 D_v: 1.740e-004
NRSolver: 9 D_u: 2.041e-001 D_v: 4.186e-008
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 10 D_u: 1.623e-001 D_v: 1.310e-007
NRSolver: 11 D_u: 1.078e-001 D_v: 2.889e-007
NRSolver: 12 D_u: 1.715e-001 D_v: 1.281e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 13 D_u: 1.997e-001 D_v: 1.462e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 14 D_u: 1.903e-001 D_v: 1.851e-007
NRSolver: 15 D_u: 1.053e-001 D_v: 4.087e-007
NRSolver: 16 D_u: 2.007e-001 D_v: 1.258e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 17 D_u: 2.114e-001 D_v: 1.792e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 18 D_u: 2.114e-001 D_v: 2.194e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 19 D_u: 2.171e-001 D_v: 2.403e-007
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 20 D_u: 2.199e-001 D_v: 2.588e-007
EngngModel info: user time consumed by solution step 20: 0.00s
Press ENTER to continue...I added X.printYouself() before the calculation of residual. My input file is as follows, now I use non-symmetric skyline:
first_contact.out
Example 16.2 from Konuyhov book
staticstructural constrainednrminiter 5 initialguess 1 manrmsteps 1 maxiter 20 nmodules 1 nsteps 20 smtype 1 solvertype 0 deltat 1 rtolf 0.001
vtkxml dofman_all tstep_all primvars 1 1
domain planestrain
outputmanager dofman_all element_all tstep_all
nbc 4 ncrosssect 1 ndofman 8 nelem 2 nic 0 nltf 3 nmat 1 nset 3
node 1 coords 2 0 0.01
node 2 coords 2 1 0.01
node 3 coords 2 1 1 bc 2 1 2
node 4 coords 2 0 1 bc 2 1 2
node 5 coords 2 -0.5 -1 bc 2 1 1
node 6 coords 2 1.5 -1 bc 2 1 1
node 7 coords 2 1.5 0
node 8 coords 2 -0.5 0
quad1planestrain 1 crosssect 1 mat 1 nlgeo 1 nodes 4 1 2 3 4
quad1planestrain 2 crosssect 1 mat 1 nlgeo 1 nodes 4 5 6 7 8
simplecs 1 thick 1
isole 1 d 1 e 1e+006 n 0.3 talpha 1
boundarycondition 1 loadtimefunction 1 prescribedvalue 0
boundarycondition 2 isimposedtimefunction 3 loadtimefunction 2 prescribedvalue -0.2
nodetosegmentfriction 3 loadtimefunction 1 masternodes 3 slavenodes 1 frictioncoefficient_d 0.1 frictioncoefficient_s 0.1 penalty_N 1e+007 penalty_T 1e+007
nodetosegmentfriction 4 loadtimefunction 1 masternodes 3 slavenodes 2 frictioncoefficient_d 0.1 frictioncoefficient_s 0.1 penalty_N 1e+007 penalty_T 1e+007
constantfunction 1 f(t) 1
piecewiselinfunction 2 f(t) 20 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 t 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
piecewiselinfunction 3 f(t) 20 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 t 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 21
set 1 nodes 1 1
set 2 nodes 1 2
set 3 nodes 2 7 8And the source file of the element in the attachment.
Do you have an idea what could be a reason for the non-convergence? The results look like the correct, but the solution does not converge.
Thanks in advance for your help.
Regards,
Vitalii
First, please correct the coding style. It takes additional mental effort for me to adapt to different coding styles all the time.
No tabs, 4 space indentation, don't indent on namespace, camel-case the functions, only classnames start with capital letters, use return values when suitable (at least for scalars)
Then, all the loops have to go. Use the functions in floatarray/floatmatrix, like add, dotProduct, beProductOf, beTProductOf, beProductTOf, plusProduct.
I will not look at any code again if you don't do this. In fact, I spent a few minutes fixing up your file, except for the loops since I don't have time for that now.
Please correct these things. I'm sure to have left some mistake in there, but this will show you how it should look (i.e. it should look like the rest of OOFEM).
The big problem, which I spotted while cleaning up the code, is due to your usage of "tStep->giveSolutionStateCounter()". This isn't what you think it is, and you should not use it.
In fact,
// It's possible to get the info about maximum substeps from the engng
history = std::vector<double>(100000, 0);
you must remove this variable. It messes up a lot of things for you. If you want the increment, or velocities, or the displacement, then you must use:
sN->giveUnknownVector(us, { D_u, D_v }, VM_Incremental, tStep, true); // or VM_Velocity for velocities
m1N->giveUnknownVector(um1, { D_u, D_v }, VM_Incremental, tStep, true);
m2N->giveUnknownVector(um2, { D_u, D_v }, VM_Incremental, tStep, true);and not try to store old solutions internally in the boundary condition. Doing that will break the simulations in the most obscure and hard-to-find ways.
But the big problem i can spot,
First, please correct the coding style. It takes additional mental effort for me to adapt to different coding styles all the time.
No tabs, 4 space indentation, don't indent on namespace, camel-case the functions, only classnames start with capital letters, use return values when suitable (at least for scalars)Then, all the loops have to go. Use the functions in floatarray/floatmatrix, like add, dotProduct, beProductOf, beTProductOf, beProductTOf, plusProduct.
I will not look at any code again if you don't do this. In fact, I spent a few minutes fixing up your file, except for the loops since I don't have time for that now.
Please correct these things. I'm sure to have left some mistake in there, but this will show you how it should look (i.e. it should look like the rest of OOFEM).The big problem, which I spotted while cleaning up the code, is due to your usage of "tStep->giveSolutionStateCounter()". This isn't what you think it is, and you should not use it.
In fact,
// It's possible to get the info about maximum substeps from the engng
history = std::vector<double>(100000, 0);
you must remove this variable. It messes up a lot of things for you. If you want the increment, or velocities, or the displacement, then you must use:sN->giveUnknownVector(us, { D_u, D_v }, VM_Incremental, tStep, true); // or VM_Velocity for velocities m1N->giveUnknownVector(um1, { D_u, D_v }, VM_Incremental, tStep, true); m2N->giveUnknownVector(um2, { D_u, D_v }, VM_Incremental, tStep, true);and not try to store old solutions internally in the boundary condition. Doing that will break the simulations in the most obscure and hard-to-find ways.
But the big problem i can spot,
Hello,
Thanks a lot for your remark. I changed the history variable and it does work like a charm.
I'll work on the programming style next week.
Regards,
Vitalii
Hello Mikael,
First, I tried to change the source code according to the OOFEM rules. You can find the outcome in the attachment, if you have any further comments let me know.
I cannot get velocities from the nodes, so far I'm just using:
double NTSF::computeConvectiveCoordinateIncrement(double ksi1, FloatMatrix xupd, Node* sN, Node *m1N,Node *m2N, TimeStep *tStep)
{
FloatArray us, um1, um2;
sN->giveUnknownVector(us, { D_u, D_v }, VM_Incremental, tStep, true); // or VM_Velocity for velocities
m1N->giveUnknownVector(um1, { D_u, D_v }, VM_Incremental, tStep, true);
m2N->giveUnknownVector(um2, { D_u, D_v }, VM_Incremental, tStep, true);
FloatMatrix xupdBefore(2, 3);
xupdBefore.at(1, 1) = xupd.at(1, 1) - um2.at(1);
xupdBefore.at(2, 1) = xupd.at(2, 1) - um2.at(2);
xupdBefore.at(1, 2) = xupd.at(1, 2) - um1.at(1);
xupdBefore.at(2, 2) = xupd.at(2, 2) - um1.at(2);
xupdBefore.at(1, 3) = xupd.at(1, 3) - us.at(1);
xupdBefore.at(2, 3) = xupd.at(2, 3) - us.at(2);
return ksi1 - computeKsi1(xupdBefore);
}Second, I still experience problems with convergence when I'm doing real calculation with a lot of nodes and significant sliding. There is no such problem in the problems without friction. I get the following outcome:
Total number of solution steps 31
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 1, metastep 1, (neq = 183)
NRSolver: Iteration ForceError DisplError
----------------------------------------------------------------------------
NRSolver: 0 D_u: *0.000e+000 0.000e+000 D_v: *0.000e+000 0.000e+000
EngngModel info: user time consumed by solution step 1: 1.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 2, metastep 1, (neq = 183)
NRSolver: Iteration ForceError DisplError
----------------------------------------------------------------------------
NRSolver: 0 D_u: *0.000e+000 0.000e+000 D_v: *1.053e+004 0.000e+000
NRSolver: 1 D_u: 2.141e-001 1.000e+000 D_v: 5.193e-001 1.000e+000
NRSolver: 2 D_u: 2.931e-001 1.140e+001 D_v: 9.575e+000 2.498e-001
NRSolver: 3 D_u: 1.915e-001 8.676e+000 D_v: 2.198e+002 3.264e+001
NRSolver: 4 D_u: 5.802e-002 1.062e+000 D_v: 5.632e-001 1.058e+000
NRSolver: 5 D_u: 2.485e-001 3.147e-002 D_v: 1.219e+002 6.217e-001
NRSolver: 6 D_u: 2.035e-001 6.577e-004 D_v: 5.782e+003 9.616e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 7 D_u: 4.560e-002 1.004e+000 D_v: 3.229e+000 1.004e+000
NRSolver: 8 D_u: 3.974e-002 1.451e+000 D_v: 1.746e+000 6.712e-001
NRSolver: 9 D_u: 2.167e-001 9.955e-001 D_v: 4.282e-001 4.911e-001
NRSolver: 10 D_u: 3.066e-001 8.796e+000 D_v: 4.453e+000 1.621e+000
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 11 D_u: 1.480e-001 2.321e-001 D_v: 4.406e+000 2.250e-001
NRSolver: 12 D_u: 2.949e-004 1.945e-001 D_v: 1.217e-001 1.825e-001
NRSolver: 13 D_u: 2.584e-005 1.125e-002 D_v: 2.659e-002 3.231e-003
NRSolver: 14 D_u: 5.891e-006 2.332e-003 D_v: 5.711e-003 7.488e-004
NRSolver: 15 D_u: 1.259e-006 4.984e-004 D_v: 1.220e-003 1.602e-004
EngngModel info: user time consumed by solution step 2: 5.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 3, metastep 1, (neq = 183)
NRSolver: Iteration ForceError DisplError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 7.148e-003 0.000e+000 D_v: 5.060e+001 0.000e+000
NRSolver: 1 D_u: 2.080e-001 9.173e-001 D_v: 5.270e-001 5.400e-001
NRSolver: 2 D_u: 3.586e-001 7.121e+000 D_v: 1.545e+001 3.037e-001
NRSolver: 3 D_u: 1.575e-001 6.580e-001 D_v: 1.733e+001 1.071e+000
NRSolver: 4 D_u: 2.254e-002 5.494e-001 D_v: 5.443e-001 2.180e-001
NRSolver: 5 D_u: 1.175e+000 3.590e+001 D_v: 4.536e+001 5.610e-001
NRSolver: 6 D_u: 1.692e-001 9.937e-001 D_v: 8.154e+000 1.105e+000
NRSolver: 7 D_u: 3.298e-002 6.398e-001 D_v: 5.422e-001 1.277e-001
NRSolver: 8 D_u: 1.152e+000 3.339e+001 D_v: 4.464e+001 5.685e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 9 D_u: 1.068e-001 9.949e-001 D_v: 2.422e+001 3.861e-001
NRSolver: 10 D_u: 2.085e-001 8.383e-001 D_v: 4.595e-001 2.415e-001
NRSolver: 11 D_u: 3.016e-001 6.186e+000 D_v: 9.415e+000 5.288e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 12 D_u: 7.351e-002 7.820e-001 D_v: 2.564e+000 1.376e-001
NRSolver: 13 D_u: 3.626e-002 8.922e-001 D_v: 3.725e-001 2.280e-001
NRSolver: 14 D_u: 8.403e-001 2.981e+001 D_v: 3.146e+001 5.204e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 15 D_u: 5.696e-002 9.826e-001 D_v: 9.482e+000 2.389e-001
NRSolver: 16 D_u: 1.394e-001 7.274e-001 D_v: 2.841e+001 1.302e+000
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 17 D_u: 1.216e-001 7.721e-001 D_v: 5.433e+000 4.023e-001
NRSolver: 18 D_u: 2.039e-001 8.412e-001 D_v: 5.599e-001 2.468e-001
NRSolver: 19 D_u: 4.272e-001 7.527e+000 D_v: 2.285e+001 3.759e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 20 D_u: 9.995e-002 6.024e-001 D_v: 1.936e+001 3.710e-001
NRSolver: 21 D_u: 2.085e-001 8.486e-001 D_v: 4.213e-001 2.469e-001
NRSolver: 22 D_u: 3.243e-001 6.187e+000 D_v: 2.544e+001 5.230e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 23 D_u: 1.012e-001 4.668e-001 D_v: 2.157e+001 5.422e-001
NRSolver: 24 D_u: 1.470e-001 1.270e-001 D_v: 2.438e+001 1.556e+000
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 25 D_u: 1.050e-001 7.801e-001 D_v: 5.253e+000 3.701e-001
NRSolver: 26 D_u: 2.033e-001 8.215e-001 D_v: 5.582e-001 2.304e-001
NRSolver: 27 D_u: 4.306e-001 7.577e+000 D_v: 2.286e+001 3.767e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 28 D_u: 1.003e-001 6.050e-001 D_v: 1.934e+001 3.710e-001
NRSolver: 29 D_u: 2.086e-001 8.485e-001 D_v: 4.215e-001 2.467e-001
NRSolver: 30 D_u: 3.055e-001 6.187e+000 D_v: 1.138e+001 2.398e-001
EngngModel info: user time consumed by solution step 3: 9.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 4, metastep 1, (neq = 183)
NRSolver: Iteration ForceError DisplError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 2.903e-001 0.000e+000 D_v: 1.546e+001 0.000e+000
NRSolver: 1 D_u: 2.204e-001 9.524e-001 D_v: 6.099e-001 3.332e-001
NRSolver: 2 D_u: 2.248e-001 4.666e+000 D_v: 9.475e+000 2.144e-001
NRSolver: 3 D_u: 1.430e-001 5.862e-001 D_v: 8.364e+001 2.808e+000
NRSolver: 4 D_u: 2.548e-001 5.200e+000 D_v: 7.719e+000 7.689e-001
NRSolver: 5 D_u: 2.231e-001 9.870e-001 D_v: 7.613e-001 3.781e-001
NRSolver: 6 D_u: 1.873e-001 4.090e+000 D_v: 7.713e+000 4.638e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 7 D_u: 1.552e-001 8.492e-001 D_v: 2.001e+000 1.786e-001
NRSolver: 8 D_u: 9.058e-002 1.553e+000 D_v: 5.167e-001 3.131e-001
NRSolver: 9 D_u: 4.183e-001 1.660e+001 D_v: 2.351e+001 5.870e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 10 D_u: 9.618e-002 1.026e+000 D_v: 1.923e+001 5.063e-001
NRSolver: 11 D_u: 2.163e-001 8.697e-001 D_v: 6.127e-001 2.653e-001
NRSolver: 12 D_u: 2.276e-001 4.735e+000 D_v: 1.109e+001 2.459e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 13 D_u: 1.240e-001 2.231e-001 D_v: 1.697e+001 5.849e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 14 D_u: 1.410e-001 8.281e-001 D_v: 2.418e+000 1.972e-001
NRSolver: 15 D_u: 7.084e-002 1.551e+000 D_v: 5.266e-001 3.122e-001
NRSolver: 16 D_u: 4.231e-001 1.710e+001 D_v: 2.356e+001 5.830e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 17 D_u: 9.669e-002 1.021e+000 D_v: 1.930e+001 5.066e-001
NRSolver: 18 D_u: 2.163e-001 8.702e-001 D_v: 6.125e-001 2.652e-001
NRSolver: 19 D_u: 2.277e-001 4.735e+000 D_v: 1.107e+001 2.455e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 20 D_u: 1.241e-001 2.233e-001 D_v: 1.694e+001 5.849e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 21 D_u: 1.411e-001 8.281e-001 D_v: 2.418e+000 1.973e-001
NRSolver: 22 D_u: 7.085e-002 1.551e+000 D_v: 5.266e-001 3.121e-001
NRSolver: 23 D_u: 4.231e-001 1.710e+001 D_v: 2.356e+001 5.831e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 24 D_u: 9.669e-002 1.021e+000 D_v: 1.930e+001 5.066e-001
NRSolver: 25 D_u: 2.163e-001 8.702e-001 D_v: 6.125e-001 2.652e-001
NRSolver: 26 D_u: 2.277e-001 4.735e+000 D_v: 1.107e+001 2.455e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 27 D_u: 1.241e-001 2.233e-001 D_v: 1.694e+001 5.849e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 28 D_u: 1.411e-001 8.281e-001 D_v: 2.418e+000 1.973e-001
NRSolver: 29 D_u: 7.085e-002 1.551e+000 D_v: 5.266e-001 3.121e-001
NRSolver: 30 D_u: 4.231e-001 1.710e+001 D_v: 2.356e+001 5.830e-001
EngngModel info: user time consumed by solution step 4: 109.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess foundDo you have a clue about the possible problem? It seems that there is fluctuations in the displacements. The solver parameters are:
staticstructural constrainednrminiter 5 initialguess 1 manrmsteps 1 maxiter 30 nmodules 1 nsteps 31 solvertype 0 deltat 1 rtolv 0.005Thanks in advance for you response.
Kind regards,
Vitalii
Hi,
it doesn't seem that you are using nonsymetric matrix storage. Your matrix become nonsymmetric with the friction so it can be the problem. For example, you can add smtype 1 to get nonsymmetric skyline.
M.
Don't use rtolv, use the rtolf parameter. Won't really make any difference, but just for future reference. rtold values are only relevant if you really know that you need them to check convergence. Else the rtolf should be fine.
I have comments on the code, but I have no suggestions for convergence problems outside of what Martin pointed out.
1. Remove _IFT_NTSF_useTangent
2. Remove GeneralBoundaryCondition::initializeFrom(ir);
3. All the code where you set "IntArray r_loc = 0" is strange. It should just be "IntArray r_loc;".
Also, other places where you set IntArrays to zero "temp_c_loc = 0" is not necessary anywhere as far as I can see. If you really needed to clear them, then "temp_c_loc.clear()" (but i don't think you actually need to clear them, as you will overwrite them anyway, or you have just created them)
4. I'd like you to rename the class to a full name, and not just an abbreviation "NTFS". Maybe "NodeEdgePenaltyContact"
5. Classes like arrays and matrices should be passed as "const FloatArray &x" if you don't plan on changing them inside the function. Else you are forcing an unnecessary copy to be made.
6. Doxygen documentation in the header file (I don't know if you already have done this)
Big things, which we don't have to have perfect before we can start including code:
1. The edges should come from an "element boundary"-set. This set contains what elements, and what boundary of what element, should be in contact. The code should then use the interpolator of this element to decide the basis functions in the boundary index given by the set.
This is what surfacetensionbc does.
The code should first be generalized to handle (almost) any basis function. Requires a bit of work, but isn't really that difficult. Just got to split the node from the rest of the N-matrix, and compound these later.
2. It should look up where the possible contacts are from a large set of elements instead of having it predefined, to allow the node(s) to slide over many edges.
These two problems we can try to work on later, and is by no means prerequisites for code inclusion.
Don't use rtolv, use the rtolf parameter. Won't really make any difference, but just for future reference. rtold values are only relevant if you really know that you need them to check convergence. Else the rtolf should be fine.
I have comments on the code, but I have no suggestions for convergence problems outside of what Martin pointed out.1. Remove _IFT_NTSF_useTangent
2. Remove GeneralBoundaryCondition::initializeFrom(ir);
3. All the code where you set "IntArray r_loc = 0" is strange. It should just be "IntArray r_loc;".
Also, other places where you set IntArrays to zero "temp_c_loc = 0" is not necessary anywhere as far as I can see. If you really needed to clear them, then "temp_c_loc.clear()" (but i don't think you actually need to clear them, as you will overwrite them anyway, or you have just created them)4. I'd like you to rename the class to a full name, and not just an abbreviation "NTFS". Maybe "NodeEdgePenaltyContact"
5. Classes like arrays and matrices should be passed as "const FloatArray &x" if you don't plan on changing them inside the function. Else you are forcing an unnecessary copy to be made.
6. Doxygen documentation in the header file (I don't know if you already have done this)
Big things, which we don't have to have perfect before we can start including code:
1. The edges should come from an "element boundary"-set. This set contains what elements, and what boundary of what element, should be in contact. The code should then use the interpolator of this element to decide the basis functions in the boundary index given by the set.
This is what surfacetensionbc does.The code should first be generalized to handle (almost) any basis function. Requires a bit of work, but isn't really that difficult. Just got to split the node from the rest of the N-matrix, and compound these later.
2. It should look up where the possible contacts are from a large set of elements instead of having it predefined, to allow the node(s) to slide over many edges.
These two problems we can try to work on later, and is by no means prerequisites for code inclusion.
Hi, Mikael.
I've fixed all "short term" issues. The cpp and header files of the element are in the attachment.
Regards,
Vitalii
Hi,
it doesn't seem that you are using nonsymetric matrix storage. Your matrix become nonsymmetric with the friction so it can be the problem. For example, you can add smtype 1 to get nonsymmetric skyline.
M.
Thanks, Martin. Indeed, I forgot to change smtype to unsymmetric skyline.
I have one more question about large strain material models. My plans are to use tutorial material (Isotropic hardening with von Mises yield locus), can I use lsmastermat to be able to work with large deformations?
You can find an example in the attachment. The intention is to work with quad1planestrain elements, instead of libeam3d.
Is it working now??
You can use lsmastermat, it just assumes stress-strain law in form S^(m) = D:(E^(m) - Ep^(m)) where
E^(m) = 1/(2m)*(U^(2m)- I) and S^(m) is work conjugated stress. Otherwise small-strain form of constitutive law is used.
So for m = 1 the constitutive law link the Green-Lagrangean stran and the second Piola-Kirchhoff stress. More reasonable, for plasticity, would by to use m = 0, i.e. the logarithmic strain and logarithmic stress.
However, I fixed some bugs in lsmastermat and it is not yet in the git version.
Is it working now??
You can use lsmastermat, it just assumes stress-strain law in form S^(m) = D:(E^(m) - Ep^(m)) where
E^(m) = 1/(2m)*(U^(2m)- I) and S^(m) is work conjugated stress. Otherwise small-strain form of constitutive law is used.So for m = 1 the constitutive law link the Green-Lagrangean stran and the second Piola-Kirchhoff stress. More reasonable, for plasticity, would by to use m = 0, i.e. the logarithmic strain and logarithmic stress.
However, I fixed some bugs in lsmastermat and it is not yet in the git version.
It calculates, and I'm going to do all verification after New Year. Can you attach file with the fixed version of lsmastermat? I'll use m = 0 for my further calculations.
If you have large rotations, but not really large strains, then you should be able to run all material models without "lsmastermat" as well.
By default, "small strain" is interpreted as the Green-Lagrange strain and the stress as 2nd Piola-Kirchoff stress.
So, for small enough strains, it'll all be the same anyway.
Yes, Mikael is right. The "lsmastermat" with m = 1 is exactly equivalent to the default large strain formulation.
I will try to make a commit of the model on Monday, or I will post the file.
The implementation is more less based on this paper
http://citeseerx.ist.psu.edu/viewdoc/do … p;type=pdf
or you can read some criticism of the model here
http://link.springer.com/chapter/10.100 … -0354-7_10
so you can decide if it is appropriate for your application.
Also, MisesMat is implemented for large deformations and is the same type of model.
Either of these options would lead to more or less the same type of material response.
Hello,
I'm working within beam01 project - the sample project. So I have two simple questions.
1. How I can add the kind of "utilities.cpp/utilities.h" to the beam01 project? I can add file to sm or core project, there you just need to add names of files in CMakeFiles.txt, but in beam01 the file is not so easy to read. Can you give me a hit how to add the file into this project? I don't have any experience with CMake, and the files are dissapearing after realoading.
2. I need information about two last steps of my calculation (node coordinates and stresses). Can I read access the data from beam01 project directly? I'm doing it so far within engngm.c, which I believe not correct way. Is there a way to read data from vtkxml?
Thanks in advance for your help.
Vitalii
1.
You add the source file to the list of files to be compiled.
Rerun cmake configuration step.
Include the header whereever you need it.
The beam01 project binary is linked against the whole of OOFEM (except main.C) so there should be no need to change that.
2.
someElement->giveIPValue(...)
someNode->giveUnknownVector(displacement, {D_u, D_v, D_w}, VM_Total, tStep)
Hi Mikael,
Thanks for your last answer, it was what I needed.
Now, all necessary components for my task have been implemented. The solution does work, however it is not really stable. The video of the example you can find here. It does work up to certain point and then it fly away, the point when it brakes depends on the mesh, but usually it is when there are more then one contact points (from the top tool) penetrates the plate. Reduction of the step reduces the problem but doesn't solve it.
Would be so kind to give hints what should I do in order to solve the problem? Is it problem due to mesh or solver?
The solver parameters are:
staticstructural constrainednrminiter 5 initialguess 1 manrmsteps 1 maxiter 30 nmodules 1 nsteps 201 smtype 1 solvertype 0 deltat 1 rtolf 0.005 The output the step before the problem and after are:
StaticStructural :: solveYourselfAt - Solving step 150, metastep 1, (neq = 363)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 9.821e-002 D_v: 3.073e+000
NRSolver: 1 D_u: 3.351e-003 D_v: 4.477e-002
NRSolver: 2 D_u: 9.360e-004 D_v: 3.634e-003
EngngModel info: user time consumed by solution step 150: 4.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 151, metastep 1, (neq = 363)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 1.021e-001 D_v: 3.032e+000
NRSolver: 1 D_u: 1.861e-001 D_v: 4.571e-001
NRSolver: 2 D_u: 8.488e-002 D_v: 2.696e-001
NRSolver: 3 D_u: 4.030e-003 D_v: 3.060e-002
NRSolver: 4 D_u: 1.079e-003 D_v: 1.005e-002
NRSolver: 5 D_u: 9.599e-005 D_v: 8.576e-004
EngngModel info: user time consumed by solution step 151: 5.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 152, metastep 1, (neq = 363)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 1.052e-001 D_v: 3.028e+000
NRSolver: 1 D_u: 1.399e-001 D_v: 4.038e-001
NRSolver: 2 D_u: 5.993e-002 D_v: 1.646e-001
NRSolver: 3 D_u: 9.194e-003 D_v: 5.531e-002
NRSolver: 4 D_u: 8.096e-003 D_v: 1.266e-002
NRSolver: 5 D_u: 1.406e-002 D_v: 6.674e-002
NRSolver: 6 D_u: 1.132e-002 D_v: 8.082e-002
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 7 D_u: 7.135e-003 D_v: 3.386e-002
NRSolver: 8 D_u: 7.153e-003 D_v: 8.888e-003
NRSolver: 9 D_u: 1.406e-002 D_v: 6.684e-002
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 10 D_u: 6.905e-003 D_v: 3.182e-002
NRSolver: 11 D_u: 7.096e-003 D_v: 8.805e-003
NRSolver: 12 D_u: 1.406e-002 D_v: 6.684e-002
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 13 D_u: 6.905e-003 D_v: 3.182e-002
NRSolver: 14 D_u: 7.096e-003 D_v: 8.805e-003
NRSolver: 15 D_u: 1.406e-002 D_v: 6.684e-002
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 16 D_u: 6.905e-003 D_v: 3.182e-002
NRSolver: 17 D_u: 7.096e-003 D_v: 8.805e-003
NRSolver: 18 D_u: 1.406e-002 D_v: 6.684e-002
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 19 D_u: 6.905e-003 D_v: 3.182e-002
NRSolver: 20 D_u: 7.096e-003 D_v: 8.805e-003
NRSolver: 21 D_u: 1.406e-002 D_v: 6.684e-002
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 22 D_u: 6.905e-003 D_v: 3.182e-002
NRSolver: 23 D_u: 7.096e-003 D_v: 8.805e-003
NRSolver: 24 D_u: 1.406e-002 D_v: 6.684e-002
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 25 D_u: 6.905e-003 D_v: 3.182e-002
NRSolver: 26 D_u: 7.096e-003 D_v: 8.805e-003
NRSolver: 27 D_u: 1.406e-002 D_v: 6.684e-002
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 28 D_u: 6.905e-003 D_v: 3.182e-002
NRSolver: 29 D_u: 7.096e-003 D_v: 8.805e-003
NRSolver: 30 D_u: 1.406e-002 D_v: 6.684e-002
EngngModel info: user time consumed by solution step 152: 30.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 153, metastep 1, (neq = 363)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 1.289e-001 D_v: 2.822e+000
NRSolver: 1 D_u: 1.183e-001 D_v: 3.955e-001
NRSolver: 2 D_u: 6.003e-002 D_v: 1.717e-001
NRSolver: 3 D_u: 1.370e-002 D_v: 1.707e-002
NRSolver: 4 D_u: 3.437e-003 D_v: 4.659e-002
NRSolver: 5 D_u: 6.615e-004 D_v: 6.193e-003
NRSolver: 6 D_u: 9.663e-005 D_v: 9.513e-004
EngngModel info: user time consumed by solution step 153: 4.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 154, metastep 1, (neq = 363)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 9.894e-002 D_v: 3.042e+000
NRSolver: 1 D_u: 5.203e-001 D_v: 6.443e-001
NRSolver: 2 D_u: 5.317e-001 D_v: 6.974e-001
NRSolver: 3 D_u: 5.338e-001 D_v: 9.929e-001
NRSolver: 4 D_u: 7.211e-001 D_v: 7.628e-001
NRSolver: 5 D_u: 6.733e-001 D_v: 9.991e-001
NRSolver: 6 D_u: 9.605e-001 D_v: 1.077e+000
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 7 D_u: 1.137e+000 D_v: 1.026e+000
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 8 D_u: 8.347e-001 D_v: 1.001e+000
NRSolver: 9 D_u: 5.514e-001 D_v: 1.251e+000
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 10 D_u: 7.661e-001 D_v: 1.114e+000
NRSolver: 11 D_u: 6.667e-001 D_v: 5.890e-001
NRSolver: 12 D_u: 1.238e+000 D_v: 1.332e+000
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 13 D_u: 9.482e-001 D_v: 1.231e+000
NRSolver: 14 D_u: 5.522e-001 D_v: 9.895e-001
NRSolver: 15 D_u: 8.657e-001 D_v: 7.892e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 16 D_u: 6.271e-001 D_v: 1.372e+000
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 17 D_u: 8.023e-001 D_v: 1.341e+000
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 18 D_u: 4.156e-001 D_v: 1.391e+000
NRSolver: 19 D_u: 4.128e-001 D_v: 7.571e-001
NRSolver: 20 D_u: 4.335e-001 D_v: 6.205e-001
NRSolver: 21 D_u: 4.455e-001 D_v: 5.687e-001
NRSolver: 22 D_u: 5.289e-001 D_v: 5.364e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 23 D_u: 3.201e-001 D_v: 4.892e-001
NRSolver: 24 D_u: 1.016e-001 D_v: 7.744e-002
NRSolver: 25 D_u: 1.144e-002 D_v: 2.412e-002
NRSolver: 26 D_u: 6.441e-004 D_v: 1.828e-003
EngngModel info: user time consumed by solution step 154: 17.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess found
StaticStructural :: solveYourselfAt - Solving step 155, metastep 1, (neq = 363)
NRSolver: Iteration ForceError
----------------------------------------------------------------------------
NRSolver: 0 D_u: 1.003e-001 D_v: 3.020e+000
NRSolver: 1 D_u: 5.740e-001 D_v: 4.011e-001
NRSolver: 2 D_u: 8.355e-001 D_v: 3.574e+000
NRSolver: 3 D_u: 8.634e-001 D_v: 6.445e-001
NRSolver: 4 D_u: 1.039e+000 D_v: 5.672e-001
NRSolver: 5 D_u: 1.252e+001 D_v: 2.466e+000
NRSolver: 6 D_u: 1.217e+000 D_v: 5.279e-001
NRSolver: 7 D_u: 1.144e+000 D_v: 5.406e-001
NRSolver: 8 D_u: 1.151e+000 D_v: 1.259e+000
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 9 D_u: 1.699e+000 D_v: 1.572e+000
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 10 D_u: 8.228e-001 D_v: 5.269e-001
NRSolver: 11 D_u: 8.339e-001 D_v: 5.312e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 12 D_u: 8.334e-001 D_v: 5.311e-001
NRSolver: 13 D_u: 8.357e-001 D_v: 5.320e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 14 D_u: 8.340e-001 D_v: 5.399e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 15 D_u: 8.352e-001 D_v: 5.401e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 16 D_u: 8.369e-001 D_v: 5.412e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 17 D_u: 8.525e-001 D_v: 6.236e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 18 D_u: 8.509e-001 D_v: 6.243e-001
NRSolver: 19 D_u: 1.033e+000 D_v: 6.602e-001
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 20 D_u: 8.413e-001 D_v: 7.866e-001
NRSolver: 21 D_u: 1.388e+000 D_v: 1.384e+000
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 22 D_u: 1.388e+000 D_v: 1.382e+000
NRSolver: 23 D_u: 1.371e+000 D_v: 1.211e+000
NRSolver: 24 D_u: 1.370e+000 D_v: 1.363e+000
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 25 D_u: 1.370e+000 D_v: 1.363e+000
Constraining increment to be 5.000000e-001 times full increment...
NRSolver: 26 D_u: 1.370e+000 D_v: 1.363e+000
NRSolver: 27 D_u: 1.347e+000 D_v: 1.301e+000
NRSolver: 28 D_u: 1.029e+000 D_v: 1.029e+000
NRSolver: 29 D_u: 1.029e+000 D_v: 1.029e+000
NRSolver: 30 D_u: 1.029e+000 D_v: 1.029e+000
EngngModel info: user time consumed by solution step 155: 25.00s
Computing initial guess
Computing old tangent
Solving for increment
Initial guess foundRegards,
Vitalii
OOFEM forum → Developers Discussion → Contact element implementation
Powered by PunBB, supported by Informer Technologies, Inc.
Currently installed 4 official extensions. Copyright © 2003–2009 PunBB.