Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Evaluation failed #985

Open
SebastienZh opened this issue May 9, 2023 · 1 comment
Open

Evaluation failed #985

SebastienZh opened this issue May 9, 2023 · 1 comment

Comments

@SebastienZh
Copy link

Hi sfepy development team!
I try to probe the cauchy stress of my FEM model. But errors occurred when I defining the evaluation of the problem.

Traceback (most recent call last):
File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\base\base.py", line 556, in getitem
ii = self.names.index(ii)
ValueError: 'ii' is not in list

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\discrete\integrals.py", line 53, in get
obj = self[name]
File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\base\base.py", line 563, in getitem
raise IndexError(msg)
IndexError: 'ii' is not in list

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "D:\PycharmWorkSpace\Box girder\FEM.py", line 129, in
stress = pb.evaluate('ev_cauchy_stress.ii.girder(m.D, u)', mode='evaluate',
File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\discrete\problem.py", line 1829, in evaluate
aux = self.create_evaluable(expression,
File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\discrete\problem.py", line 1771, in create_evaluable
out = create_evaluable(expression, self.fields, materials,
File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\discrete\evaluate.py", line 230, in create_evaluable
equations = Equations.from_conf({'tmp' : expression},
File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\discrete\equations.py", line 65, in from_conf
eq = Equation.from_desc(name, desc, variables, regions,
File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\discrete\equations.py", line 743, in from_desc
terms = Terms.from_desc(term_descs, regions, integrals)
File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\terms\terms.py", line 218, in from_desc
term = Term.from_desc(constructor, td, region, integrals=integrals)
File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\terms\terms.py", line 354, in from_desc
integral = integrals.get(desc.integral)
File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\discrete\integrals.py", line 56, in get
raise ValueError('integral %s is not defined!' % name)
ValueError: integral ii is not defined!

It says the integral ii is not defined. But I have defined the integral earlier.
integral = Integral('ii', order=3)

The full code is listed below:

from __future__ import absolute_import
import numpy as np
import pandas as pd
from sfepy.base.base import IndexedStruct
from sfepy.discrete import (FieldVariable, Material, Integral, Function, Equation, Equations, Problem)
from sfepy.discrete.conditions import Conditions, EssentialBC
from sfepy.discrete.fem import Mesh, FEDomain, Field
from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson
from sfepy.solvers.ls import ScipyDirect
from sfepy.solvers.nls import Newton
from sfepy.terms import Term
from sfepy.base.base import Struct
from sfepy.postprocess.probes_vtk import Probe
from six.moves import range

hf = 1
wf = 2
hr = 5
wr = 1
hs = 1
ws = 6
theta = 0.25 * np.pi
span_len = 2
p = -0.5 * 28.14 / (wf + (hr + hs) / np.tan(theta) + ws * 0.5)  


def get_bearing_point(coors, domain=None):
    coord_df = pd.DataFrame(coors)
    bottom_coord = coord_df[(coord_df[2] == coord_df[2].min()) & (coord_df[1] == coord_df[1].min())]  # 梁端底面
    min_x = bottom_coord[0].min()
    max_x = bottom_coord[0].max()
    x_1_4 = (max_x - min_x) / 4 + min_x
    x_3_4 = 3 * (max_x - min_x) / 4 + min_x
    bearing_point_1 = bottom_coord[abs(bottom_coord[0] - x_1_4) < 0.5]
    bearing_point_2 = bottom_coord[abs(bottom_coord[0] - x_3_4) < 0.5]
    flag = np.array(
        list(set(
            list(bearing_point_1.index) + list(bearing_point_2.index)
        ))
    )
    return flag


def linear_compress(ts, coor, pressure, mode=None, **kwargs):
    if mode == 'qp':
        val = np.tile(-1 * pressure, (coor.shape[0], 1, 1))

        return {'val': val}


pressure = Function('linear_compress', linear_compress, extra_args={'pressure': p})

functions = {
    'get_bearing_point': (get_bearing_point,)
}

mesh = Mesh.from_file('box_girder.vtk')

domain = FEDomain('domain', mesh)

bounding_box = domain.get_mesh_bounding_box()
min_y, max_y = bounding_box[:, 1]
min_z, max_z = bounding_box[:, 2]
eps_y = 1e-8 * (max_y - min_y)
eps_z = 1e-8 * (max_z - min_z)

bearing = domain.create_region('bearing', 'vertices by get_bearing_point')  # 梁端固定支座
bearing = domain.create_region('bearing', 'vertices in y < %.10f' % (min_y + eps_y), 'facet')
deck = domain.create_region('deck', 'vertices in z > %.10f' % (max_z - eps_z), 'facet')  # 梁顶面
girder = domain.create_region('girder', 'all')

field = Field.from_args('fu', np.float64, 'vector', girder, approx_order=2)
u = FieldVariable('u', 'unknown', field)
v = FieldVariable('v', 'test', field, primary_var_name='u')

m = Material('m', D=stiffness_from_youngpoisson(3, 3.5e10, 0.2))
f = Material('f', val=[[0], [0], [-23520]])
load = Material('load', function=pressure)

integral = Integral('ii', order=3)

t1 = Term.new('dw_lin_elastic(m.D, v, u)', integral, girder, m=m, v=v, u=u)  # 弹性能积分项
t2 = Term.new('dw_volume_lvf(f.val, v)', integral, girder, f=f, v=v)  # 体积力积分项
t3 = Term.new('dw_surface_ltr(load.val, v)', integral, deck, load=load, v=v)  # 面力积分项

eq = Equation('balance', t1 - t2 - t3)
eqs = Equations([eq])

fix_u = EssentialBC('fix_u', bearing, {'u.all': 0.0})

ls = ScipyDirect({})
nls_status = IndexedStruct()
nls = Newton({'i_max': 1, 'eps_a': 1e-10}, lin_solver=ls, status=nls_status)

pb = Problem('elasticity', equations=eqs)
pb.save_regions_as_groups('regions')

pb.set_bcs(ebcs=Conditions([fix_u]))
pb.set_solver(nls)

status = IndexedStruct()
vec = pb.solve(status=status)

print('Nonlinear solver status:\n', nls_status)
print('Stationary solver status:\n', status)

pb.save_state('linear_elasticity.vtk', vec)


out = {}  
stress = pb.evaluate('ev_cauchy_stress.ii.girder(m.D, u)', mode='evaluate')
out['cauchy_stress'] = Struct(name='output_data', mode='cell', data=stress, dofs=None)

I would greatly appreciate your help of any kind.
Thank you!

@SebastienZh
Copy link
Author

By the way is there any tutorial that shows how to write postprocess code interactively? I believe there is only tutorials of the problem defining part.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant