Works like a charm. Thank you again.
Roman
2015-10-06 17:49 GMT+03:00 David Knezevic <david.kneze...@akselos.com>:
> On Tue, Oct 6, 2015 at 10:46 AM, Вася Бубликов <gream...@gmail.com> wrote:
>
>> David
>>
>>
>> It pains me to trouble you again but
>> dynamic_cast<PetscMatrix<Number>*>(system.matrix);
>> returns NULL!
>>
>
>
> Make sure you do the cast after equation_systems.init();
>
> David
>
>
>
>
>
>> 2015-10-06 17:11 GMT+03:00 David Knezevic <david.kneze...@akselos.com>:
>>
>>> Please reply to the list as well.
>>>
>>> Instead of creating a new matrix, you should cast system.matrix to a
>>> PetscMatrix and then call MatSetOption on that.
>>>
>>> David
>>>
>>>
>>>
>>>
>>> On Tue, Oct 6, 2015 at 10:06 AM, Вася Бубликов <gream...@gmail.com>
>>> wrote:
>>>
>>>> David
>>>>
>>>> Thank you for your answer. Here is a short program to try out your
>>>> first suggestion:
>>>>
>>>> #include <iostream>
>>>> #include <vector>
>>>> #include <set>
>>>> #include <algorithm>
>>>> #include <fstream>
>>>> #include <string>
>>>> #include <ctime>
>>>> // libmesh includes
>>>> #include "libmesh/libmesh.h"
>>>> #include "libmesh/mesh.h"
>>>> #include "libmesh/linear_implicit_system.h"
>>>> #include "libmesh/equation_systems.h"
>>>> #include "libmesh/fe.h"
>>>> #include "libmesh/quadrature_trap.h"
>>>> #include "libmesh/dof_map.h"
>>>> #include "libmesh/sparse_matrix.h"
>>>> #include "libmesh/petsc_matrix.h"
>>>> #include "libmesh/numeric_vector.h"
>>>> #include "libmesh/dense_matrix.h"
>>>> #include "libmesh/dense_submatrix.h"
>>>> #include "libmesh/dense_vector.h"
>>>> #include "libmesh/dense_subvector.h"
>>>> #include "libmesh/perf_log.h"
>>>> #include "libmesh/elem.h"
>>>> #include "libmesh/boundary_info.h"
>>>> #include "libmesh/string_to_enum.h"
>>>> #include "libmesh/getpot.h"
>>>> #include "libmesh/vtk_io.h"
>>>> #include "libmesh/mesh_data.h"
>>>> #include "libmesh/mesh_base.h"
>>>>
>>>> using namespace libMesh;
>>>>
>>>> int main (int argc, char** argv)
>>>> {
>>>> LibMeshInit init (argc, argv);
>>>> SerialMesh mesh(init.comm());
>>>> mesh.read("SquareWithCircle.unv");
>>>> mesh.print_info();
>>>> EquationSystems equation_systems (mesh);
>>>> LinearImplicitSystem& system =
>>>> equation_systems.add_system<LinearImplicitSystem> ("Elasticity");
>>>> PetscMatrix<Number> Matrix(init.comm());
>>>> std::cout << "Enable non neighbours interactions" << std::endl;
>>>> MatSetOption( Matrix.mat(), MAT_NEW_NONZERO_ALLOCATION_ERR,
>>>> PETSC_FALSE);
>>>> unsigned int u_var = system.add_variable("u", FIRST, LAGRANGE); //
>>>> displacement X
>>>> unsigned int v_var = system.add_variable("v", FIRST, LAGRANGE); //
>>>> displacement Y
>>>> equation_systems.init();
>>>> equation_systems.print_info();
>>>> Matrix.attach_dof_map(system.get_dof_map());
>>>> std::cout << "Initialize matrix" << std::endl;
>>>> Matrix.init();
>>>> std::cout << "Make my matrix the system matrix" << std::endl;
>>>> system.matrix = &Matrix;
>>>> return 0;
>>>> }
>>>>
>>>>
>>>> it compiles fine but then I run it MatSetOption command causes an error:
>>>>
>>>> [0]PETSC ERROR: --------------------- Error Message
>>>> --------------------------------------------------------------
>>>> [0]PETSC ERROR: Corrupt argument:
>>>> http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
>>>> [0]PETSC ERROR: Invalid Pointer to Object: Parameter # 1
>>>> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
>>>> for trouble shooting.
>>>> [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015
>>>> [0]PETSC ERROR: ./error.o on a arch-linux-c-opt named izmailovo by
>>>> gream Tue Oct 6 17:01:25 2015
>>>> [0]PETSC ERROR: Configure options --prefix=/opt/petsc/linux-c-opt
>>>> --PETSC_ARCH=arch-linux-c-opt --with-shared-libraries=1 --COPTFLAGS=-O3
>>>> --CXXOPTFLAGS=-O3 --with-mumps=1 --with-hdf5=1 --with-scalapack=1
>>>> --with-suitesparse=1 --with-metis=1 --with-parmetis=1 --with-ptscotch=1
>>>> --with-ptscotch-lib="[libesmumps.so,libptscotch.so,libptscotcherr.so,libscotch.so,libscotcherr.so,libbz2.so]"
>>>> --with-ptscotch-include=/usr/include/scotch --with-superlu=1
>>>> --with-superlu-lib=-lsuperlu --with-superlu-include=/usr/include/superlu
>>>> --with-ml=1
>>>> [0]PETSC ERROR: #1 MatSetOption() line 5261 in
>>>> /tmp/yaourt-tmp-gream/aur-petsc/src/petsc-3.6.1/src/mat/interface/matrix.c
>>>>
>>>> Roman
>>>>
>>>> 2015-10-06 13:20 GMT+03:00 David Knezevic <david.kneze...@akselos.com>:
>>>>
>>>>> This is a PETSc error that you're hitting because you're adding new
>>>>> non-zeros in the matrix. PETSc tries to help you by throwing an error if
>>>>> you do this, because allocating new non-zeros (outside the specified
>>>>> sparsity pattern) can slow down your program a lot.
>>>>>
>>>>> The easy fix is indicated in the PETSc error message:
>>>>> "Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to
>>>>> turn
>>>>> off this check"
>>>>> To do this, you need to cast your matrix to a PetscMatrix, then do
>>>>> petsc_matrix->mat() to get the PETSc Mat that you can pass to
>>>>> MatSetOption.
>>>>>
>>>>> The more involved fix is to augment your sparsity pattern (e.g. see
>>>>> miscellaneous_ex9).
>>>>>
>>>>> David
>>>>>
>>>>>
>>>>>
>>>>> On Tue, Oct 6, 2015 at 6:13 AM, Вася Бубликов <gream...@gmail.com>
>>>>> wrote:
>>>>>
>>>>>> Dear Sirs,
>>>>>>
>>>>>> I'm trying to make Libmesh work with a simple linear elasticity
>>>>>> problem
>>>>>> using linear systems example 4 as an example. I need to force the top
>>>>>> boundary remain flat, i.e. all normal displacements are of the same
>>>>>> value
>>>>>> to be computed.
>>>>>> To obtain this, for example to make the displacements for DOFs
>>>>>> i1,i2,i3..in
>>>>>> the same as for the DOF j, I add penalty (1e10) to K(in,in) and
>>>>>> subtract it
>>>>>> from K(in,j) meaning to add 1e10*(U_in - U_j) = 0 to the in'th
>>>>>> equation
>>>>>> in the assembly function. But as soon as in and j are on different
>>>>>> finite
>>>>>> elements I get the same error no matter which particular c++ function
>>>>>> I use
>>>>>> (system.matrix->add(i,j,value)
>>>>>> system.matrix->add_matrix(DenseMatrix)
>>>>>> or system.matrix->add(SparseMatrix) ). I would be very mush obliged
>>>>>> if you
>>>>>> helped me.
>>>>>>
>>>>>> Kind regards,
>>>>>> Roman Khudobin
>>>>>> Research Assistant ICHPH RAS
>>>>>>
>>>>>> [0]PETSC ERROR: --------------------- Error Message
>>>>>> --------------------------------------------------------------
>>>>>> [0]PETSC ERROR: Argument out of range
>>>>>> [0]PETSC ERROR: New nonzero at (17,47) caused a malloc
>>>>>> Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to
>>>>>> turn
>>>>>> off this check
>>>>>> [0]PETSC ERROR: See
>>>>>> http://www.mcs.anl.gov/petsc/documentation/faq.html for
>>>>>> trouble shooting.
>>>>>> [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015
>>>>>> [0]PETSC ERROR: ./constrain.o on a arch-linux-c-opt named izmailovo by
>>>>>> gream Tue Oct 6 12:55:21 2015
>>>>>> [0]PETSC ERROR: Configure options --prefix=/opt/petsc/linux-c-opt
>>>>>> --PETSC_ARCH=arch-linux-c-opt --with-shared-libraries=1
>>>>>> --COPTFLAGS=-O3
>>>>>> --CXXOPTFLAGS=-O3 --with-mumps=1 --with-hdf5=1 --with-scalapack=1
>>>>>> --with-suitesparse=1 --with-metis=1 --with-parmetis=1
>>>>>> --with-ptscotch=1
>>>>>>
>>>>>> --with-ptscotch-lib="[libesmumps.so,libptscotch.so,libptscotcherr.so,libscotch.so,libscotcherr.so,libbz2.so]"
>>>>>> --with-ptscotch-include=/usr/include/scotch --with-superlu=1
>>>>>> --with-superlu-lib=-lsuperlu
>>>>>> --with-superlu-include=/usr/include/superlu
>>>>>> --with-ml=1
>>>>>> [0]PETSC ERROR: #1 MatSetValues_SeqAIJ() line 485 in
>>>>>>
>>>>>> /tmp/yaourt-tmp-gream/aur-petsc/src/petsc-3.6.1/src/mat/impls/aij/seq/aij.c
>>>>>> [0]PETSC ERROR: #2 MatSetValues() line 1173 in
>>>>>>
>>>>>> /tmp/yaourt-tmp-gream/aur-petsc/src/petsc-3.6.1/src/mat/interface/matrix.c
>>>>>>
>>>>>> --------------------------------------------------------------------------
>>>>>> MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
>>>>>> with errorcode 1.
>>>>>>
>>>>>> NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
>>>>>> You may or may not see output from other processes, depending on
>>>>>> exactly when Open MPI kills them.
>>>>>>
>>>>>> --------------------------------------------------------------------------
>>>>>>
>>>>>>
>>>>>> ------------------------------------------------------------------------------
>>>>>>
>>>>>> _______________________________________________
>>>>>> Libmesh-users mailing list
>>>>>> Libmesh-users@lists.sourceforge.net
>>>>>> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>>>>>
>>>>>>
>>>>>
>>>>
>>>
>>
------------------------------------------------------------------------------
_______________________________________________
Libmesh-users mailing list
Libmesh-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-users
