David

It pains me to trouble you again but
dynamic_cast<PetscMatrix<Number>*>(system.matrix);
returns NULL!

Roman

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

Reply via email to