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
