Python Problem.save_regions_as_groups方法代码示例

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import save_regions_as_groups [as 别名]
def main():
    from sfepy import data_dir

    parser = OptionParser(usage=usage, version="%prog")
    parser.add_option("-s", "--show", action="store_true", dest="show", default=False, help=help["show"])
    options, args = parser.parse_args()

    mesh = Mesh.from_file(data_dir + "/meshes/2d/rectangle_tri.mesh")
    domain = Domain("domain", mesh)

    min_x, max_x = domain.get_mesh_bounding_box()[:, 0]
    eps = 1e-8 * (max_x - min_x)
    omega = domain.create_region("Omega", "all")
    gamma1 = domain.create_region("Gamma1", "vertices in x < %.10f" % (min_x + eps), "facet")
    gamma2 = domain.create_region("Gamma2", "vertices in x > %.10f" % (max_x - eps), "facet")

    field = Field.from_args("fu", nm.float64, "vector", omega, approx_order=2)

    u = FieldVariable("u", "unknown", field)
    v = FieldVariable("v", "test", field, primary_var_name="u")

    m = Material("m", lam=1.0, mu=1.0)
    f = Material("f", val=[[0.02], [0.01]])

    integral = Integral("i", order=3)

    t1 = Term.new("dw_lin_elastic_iso(m.lam, m.mu, v, u)", integral, omega, m=m, v=v, u=u)
    t2 = Term.new("dw_volume_lvf(f.val, v)", integral, omega, f=f, v=v)
    eq = Equation("balance", t1 + t2)
    eqs = Equations([eq])

    fix_u = EssentialBC("fix_u", gamma1, {"u.all": 0.0})

    bc_fun = Function("shift_u_fun", shift_u_fun, extra_args={"shift": 0.01})
    shift_u = EssentialBC("shift_u", gamma2, {"u.0": bc_fun})

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = Problem("elasticity", equations=eqs, nls=nls, ls=ls)
    pb.save_regions_as_groups("regions")

    pb.time_update(ebcs=Conditions([fix_u, shift_u]))

    vec = pb.solve()
    print nls_status

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

    if options.show:
        view = Viewer("linear_elasticity.vtk")
        view(vector_mode="warp_norm", rel_scaling=2, is_scalar_bar=True, is_wireframe=True)

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import save_regions_as_groups [as 别名]

#.........这里部分代码省略.........
    eq1 = Equation('stiffness', t1)
    eq2 = Equation('mass', t2)
    lhs_eqs = Equations([eq1, eq2])

    pb = Problem('modal', equations=lhs_eqs)

    if options.bc_kind == 'free':
        pb.time_update()
        n_rbm = dim * (dim + 1) / 2

    elif options.bc_kind == 'cantilever':
        fixed = EssentialBC('Fixed', bottom, {'u.all' : 0.0})
        pb.time_update(ebcs=Conditions([fixed]))
        n_rbm = 0

    elif options.bc_kind == 'fixed':
        fixed = EssentialBC('Fixed', bottom_top, {'u.all' : 0.0})
        pb.time_update(ebcs=Conditions([fixed]))
        n_rbm = 0

    else:
        raise ValueError('unsupported BC kind! (%s)' % options.bc_kind)

    if options.ignore is not None:
        n_rbm = options.ignore

    pb.update_materials()

    # Assemble stiffness and mass matrices.
    mtx_k = eq1.evaluate(mode='weak', dw_mode='matrix', asm_obj=pb.mtx_a)
    mtx_m = mtx_k.copy()
    mtx_m.data[:] = 0.0
    mtx_m = eq2.evaluate(mode='weak', dw_mode='matrix', asm_obj=mtx_m)

    try:
        eigs, svecs = eig_solver(mtx_k, mtx_m, options.n_eigs + n_rbm,
                                 eigenvectors=True)

    except sla.ArpackNoConvergence as ee:
        eigs = ee.eigenvalues
        svecs = ee.eigenvectors
        output('only %d eigenvalues converged!' % len(eigs))

    output('%d eigenvalues converged (%d ignored as rigid body modes)' %
           (len(eigs), n_rbm))

    eigs = eigs[n_rbm:]
    svecs = svecs[:, n_rbm:]

    omegas = nm.sqrt(eigs)
    freqs = omegas / (2 * nm.pi)

    output('number |         eigenvalue |  angular frequency '
           '|          frequency')
    for ii, eig in enumerate(eigs):
        output('%6d | %17.12e | %17.12e | %17.12e'
               % (ii + 1, eig, omegas[ii], freqs[ii]))

    # Make full eigenvectors (add DOFs fixed by boundary conditions).
    variables = pb.get_variables()

    vecs = nm.empty((variables.di.ptr[-1], svecs.shape[1]),
                    dtype=nm.float64)
    for ii in range(svecs.shape[1]):
        vecs[:, ii] = variables.make_full_vec(svecs[:, ii])

    # Save the eigenvectors.
    out = {}
    state = pb.create_state()
    for ii in range(eigs.shape[0]):
        state.set_full(vecs[:, ii])
        aux = state.create_output_dict()
        strain = pb.evaluate('ev_cauchy_strain.i.Omega(u)',
                             integrals=Integrals([integral]),
                             mode='el_avg', verbose=False)
        out['u%03d' % ii] = aux.popitem()[1]
        out['strain%03d' % ii] = Struct(mode='cell', data=strain)

    pb.save_state('eigenshapes.vtk', out=out)
    pb.save_regions_as_groups('regions')

    if len(eigs) and options.show:
        # Show the solution. If the approximation order is greater than 1, the
        # extra DOFs are simply thrown away.
        from sfepy.postprocess.viewer import Viewer
        from sfepy.postprocess.domain_specific import DomainSpecificPlot

        scaling = 0.05 * dims.max() / nm.abs(vecs).max()

        ds = {}
        for ii in range(eigs.shape[0]):
            pd = DomainSpecificPlot('plot_displacements',
                                    ['rel_scaling=%s' % scaling,
                                     'color_kind="tensors"',
                                     'color_name="strain%03d"' % ii])
            ds['u%03d' % ii] = pd

        view = Viewer('eigenshapes.vtk')
        view(domain_specific=ds, only_names=sorted(ds.keys()),
             is_scalar_bar=False, is_wireframe=True)

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import save_regions_as_groups [as 别名]

#.........这里部分代码省略.........
    gamma5 = domain.create_region('Bottom',
                                  'vertices in z < %.10f' % (min_z + eps),
                                  'facet')
    gamma6 = domain.create_region('Top',
                                  'vertices in z > %.10f' % (max_z - eps),
                                  'facet')



    field = Field.from_args('fu', nm.float64, 'vector', omega, approx_order=2)

    u = FieldVariable('u', 'unknown', field)
    v = FieldVariable('v', 'test', field, primary_var_name='u')
    mu=1.1
    lam=1.0
    m = Material('m', lam=lam, mu=mu)
    f = Material('f', val=[[0.0], [0.0],[-1.0]])
    load = Material('Load',val=[[0.0],[0.0],[-Load]])

    D = stiffness_from_lame(3,lam, mu)
    mat = Material('Mat', D=D)

    get_mat = Function('get_mat1',get_mat1)
    get_mat_f = Function('get_mat_f',get_mat1)

    integral = Integral('i', order=3)
    s_integral = Integral('is',order=2)

    t1 = Term.new('dw_lin_elastic(Mat.D, v, u)',
         integral, omega, Mat=mat, v=v, u=u)
    t2 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, f=f, v=v)
    #t3 = Term.new('DotProductSurfaceTerm(Load.val, v)',s_integral,gamma5,Load=load,v=v)
    t3 = Term.new('dw_surface_ltr( Load.val, v )',s_integral,gamma6,Load=load,v=v)
    eq = Equation('balance', t1 + t2 + t3)
    eqs = Equations([eq])

    fix_u = EssentialBC('fix_u', gamma1, {'u.all' : 0.0})
    left_bc  = EssentialBC('Left',  gamma1, {'u.0' : 0.0})
    right_bc = EssentialBC('Right', gamma2, {'u.0' : 0.0})
    back_bc = EssentialBC('Front', gamma3, {'u.1' : 0.0})
    front_bc = EssentialBC('Back', gamma4, {'u.1' : 0.0})
    bottom_bc = EssentialBC('Bottom', gamma5, {'u.all' : 0.0})
    top_bc = EssentialBC('Top', gamma6, {'u.2' : 0.2})

    bc=[left_bc,right_bc,back_bc,front_bc,bottom_bc]
    #bc=[bottom_bc,top_bc]


    ##############################
    #  ##### SOLVER SECTION  #####
    ##############################

    conf = Struct(method='bcgsl', precond='jacobi', sub_precond=None,
                  i_max=10000, eps_a=1e-50, eps_r=1e-10, eps_d=1e4,
                  verbose=True)

    ls = PETScKrylovSolver(conf)

    file.write(str(ls.name)+' '+str(ls.conf.method)+' '+str(ls.conf.precond)+' '+str(ls.conf.eps_r)+' '+str(ls.conf.i_max)+'\n' )

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

    pb = Problem('elasticity', equations=eqs, nls=nls, ls=ls)

    dd=pb.get_materials()['Mat']
    dd.set_function(get_mat1)
    #xload = pb.get_materials()['f']
    #xload.set_function(get_mat_f)

    pb.save_regions_as_groups('regions')

    pb.time_update(ebcs=Conditions(bc))

    vec = pb.solve()
    print nls_status


    file.write('TIME TO SOLVE\n')
    file.write(str(nls.status.time_stats['solve'])+'\n')
    file.write('TIME TO CREATE MATRIX\n')
    file.write(str(nls.status.time_stats['matrix'])+'\n')

    ev = pb.evaluate
    out = vec.create_output_dict()
    strain = ev('ev_cauchy_strain.3.Omega(u)', mode='el_avg')
    stress = ev('ev_cauchy_stress.3.Omega(Mat.D, u)', mode='el_avg',
                copy_materials=False)

    out['cauchy_strain'] = Struct(name='output_data', mode='cell',
                                  data=strain, dofs=None)
    out['cauchy_stress'] = Struct(name='output_data', mode='cell',
                                  data=stress, dofs=None)

    pb.save_state('strain.vtk', out=out)

    print nls_status


    file.close()

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import save_regions_as_groups [as 别名]
def main():
    from sfepy import data_dir

    parser = ArgumentParser()
    parser.add_argument('--version', action='version', version='%(prog)s')
    parser.add_argument('-s', '--show',
                        action="store_true", dest='show',
                        default=False, help=helps['show'])
    options = parser.parse_args()

    mesh = Mesh.from_file(data_dir + '/meshes/2d/rectangle_tri.mesh')
    domain = FEDomain('domain', mesh)

    min_x, max_x = domain.get_mesh_bounding_box()[:,0]
    eps = 1e-8 * (max_x - min_x)
    omega = domain.create_region('Omega', 'all')
    gamma1 = domain.create_region('Gamma1',
                                  'vertices in x < %.10f' % (min_x + eps),
                                  'facet')
    gamma2 = domain.create_region('Gamma2',
                                  'vertices in x > %.10f' % (max_x - eps),
                                  'facet')

    field = Field.from_args('fu', nm.float64, 'vector', omega,
                            approx_order=2)

    u = FieldVariable('u', 'unknown', field)
    v = FieldVariable('v', 'test', field, primary_var_name='u')

    m = Material('m', D=stiffness_from_lame(dim=2, lam=1.0, mu=1.0))
    f = Material('f', val=[[0.02], [0.01]])

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

    t1 = Term.new('dw_lin_elastic(m.D, v, u)',
                  integral, omega, m=m, v=v, u=u)
    t2 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, f=f, v=v)
    eq = Equation('balance', t1 + t2)
    eqs = Equations([eq])

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

    bc_fun = Function('shift_u_fun', shift_u_fun,
                      extra_args={'shift' : 0.01})
    shift_u = EssentialBC('shift_u', gamma2, {'u.0' : bc_fun})

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

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

    pb.time_update(ebcs=Conditions([fix_u, shift_u]))

    vec = pb.solve()
    print(nls_status)

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

    if options.show:
        view = Viewer('linear_elasticity.vtk')
        view(vector_mode='warp_norm', rel_scaling=2,
             is_scalar_bar=True, is_wireframe=True)

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import save_regions_as_groups [as 别名]
    def _solve(self, property_array):
        """
        Solve the Sfepy problem for one sample.

        Args:
          property_array: array of shape (Nx, Ny, 2) where the last
          index is for Lame's parameter and shear modulus,
          respectively.

        Returns:
          the strain field of shape (Nx, Ny, 2) where the last
          index represents the x and y displacements

        """
        shape = property_array.shape[:-1]
        mesh = self._get_mesh(shape)
        domain = Domain('domain', mesh)

        region_all = domain.create_region('region_all', 'all')

        field = Field.from_args('fu', np.float64, 'vector', region_all,
                                approx_order=2)

        u = FieldVariable('u', 'unknown', field)
        v = FieldVariable('v', 'test', field, primary_var_name='u')

        m = self._get_material(property_array, domain)

        integral = Integral('i', order=4)

        t1 = Term.new('dw_lin_elastic_iso(m.lam, m.mu, v, u)',
                      integral, region_all, m=m, v=v, u=u)
        eq = Equation('balance_of_forces', t1)
        eqs = Equations([eq])

        epbcs, functions = self._get_periodicBCs(domain)
        ebcs = self._get_displacementBCs(domain)
        lcbcs = self._get_linear_combinationBCs(domain)

        ls = ScipyDirect({})

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

        pb.time_update(ebcs=ebcs, epbcs=epbcs, lcbcs=lcbcs, functions=functions)

        ev = pb.get_evaluator()
        nls = Newton({}, lin_solver=ls,
                     fun=ev.eval_residual, fun_grad=ev.eval_tangent_matrix)

        pb.set_solvers_instances(ls, nls)

        vec = pb.solve()

        u = vec.create_output_dict()['u'].data
        u_reshape = np.reshape(u, (tuple(x + 1 for x in shape) + u.shape[-1:]))

        dims = domain.get_mesh_bounding_box().shape[1]
        strain = np.squeeze(pb.evaluate('ev_cauchy_strain.{dim}.region_all(u)'.format(dim=dims),
                                        mode='el_avg'))
        strain_reshape = np.reshape(strain, (shape + strain.shape[-1:]))

        stress = np.squeeze(pb.evaluate('ev_cauchy_stress.{dim}.region_all(m.D, u)'.format(dim=dims),
                                        mode='el_avg'))
        stress_reshape = np.reshape(stress, (shape + stress.shape[-1:]))

        return strain_reshape, u_reshape, stress_reshape

# 需要导入模块: from sfepy.discrete import Problem [as 别名]
# 或者: from sfepy.discrete.Problem import save_regions_as_groups [as 别名]

#.........这里部分代码省略.........
#if hasattr(ls.name,'ls.scipy_iterative'):
    file.write(str(ls.name)+' '+str(ls.conf.method)+' '+str(ls.conf.precond)+' '+str(ls.conf.eps_r)+' '+str(ls.conf.i_max)+'\n' )
        #    else:
#file.write(str(ls.name)+' '+str(ls.conf.method)+'\n')



   
   
   # conf = Struct(method='bcgsl', precond='jacobi', sub_precond=None,
   #                 i_max=10000, eps_a=1e-50, eps_r=1e-8, eps_d=1e4,
#              verbose=True)
            
                 
                 
#ls = PETScKrylovSolver(conf)



#ls = ScipyIterative({'method':'bicgstab','i_max':100,'eps_r':1e-10})


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

    pb = Problem('elasticity', equations=eqs, nls=nls, ls=ls)
    


    dd=pb.get_materials()['Mat']
    dd.set_function(get_mat1)
    
    
    pb.save_regions_as_groups('regions')

    #pb.time_update(ebcs=Conditions([fix_u, shift_u]))

    pb.time_update(ebcs=Conditions(bc))
    pb.save_regions_as_groups('regions')

#ls = ScipyIterative({'method':'bicgstab','i_max':100,'eps_r':1e-10})


#   A = pb.mtx_a
#   M = spilu(A,fill_factor = 1)
    
    #conf = Struct(solvers ='ScipyIterative',method='bcgsl', sub_precond=None,
# i_max=1000, eps_r=1e-8)
        
#pb.set_conf_solvers(conf)
    vec = pb.solve()
    print nls_status
    file.write('TIME TO SOLVE\n')
    file.write(str(nls.status.time_stats['solve'])+'\n')
    file.write('TIME TO CREATE MATRIX\n')
    file.write(str(nls.status.time_stats['matrix'])+'\n')
    #out = post_process(out, pb, state, extend=False)
    ev = pb.evaluate
    out = vec.create_output_dict()
    strain = ev('ev_cauchy_strain.3.Omega(u)', mode='el_avg')
    stress = ev('ev_cauchy_stress.3.Omega(Mat.D, u)', mode='el_avg',
                copy_materials=False)

    out['cauchy_strain'] = Struct(name='output_data', mode='cell',
                                  data=strain, dofs=None)
    out['cauchy_stress'] = Struct(name='output_data', mode='cell',