Hello, GetFEM project.
I have a problem using the add_bilaplacian_brick method of the Model
object. My wish is to get the following stiffness matrix in the attached
script file, but I got an empty stiffness matrix.
Trace 2 in getfem_fourth_order.cc, line 88: Stiffness matrix assembly of a
bilaplacian term
Computed matrix is...
matrix(2, 2)
( )
( )
Expected matrix is...
matrix(2, 2)
( (r0, 100) (r1, -100) )
( (r0, -100) (r1, 100) )
I'm sorry, but could you please review my script to see if there is any
problem?
I used the developing version in the git repository.
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Python GetFEM interface
#
# Copyright (C) 2020-2020 Tetsuo Koyama.
#
# This file is a part of GetFEM
#
# GetFEM is free software; you can redistribute it and/or modify it
# under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation; either version 2.1 of the License, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
# License for more details.
# You should have received a copy of the GNU Lesser General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
#
############################################################################
""" 1D Truss problem test.
This program is used to check that python-getfem is working. This is
also a good example of use of GetFEM.
$Id$
"""
import getfem as gf
import numpy as np
E = 1.0E+05 # Young Modulus (N/mm^2)
A = 1.0 # Cross-Sectional Area (mm^2)
L = 1000.0 # Length (mm)
mesh = gf.Mesh("cartesian", [0, L])
mfu = gf.MeshFem(mesh, 1)
element_degree = 1
mfu.set_classical_fem(element_degree)
mim = gf.MeshIm(mesh, element_degree*2)
model = gf.Model("real")
model.add_fem_variable("u", mfu)
model.add_initialized_data("D", [E*A])
model.add_bilaplacian_brick(mim, "u", "D")
model.assembly()
SM = model.tangent_matrix()
print("Computed matrix is...")
print(SM)
print("Expected matrix is...")
SM = gf.Spmat("empty", 2, 2)
SM.add(0, 0, E*A/L)
SM.add(0, 1, -E*A/L)
SM.add(1, 0, -E*A/L)
SM.add(1, 1, E*A/L)
print(SM)
