From: Rafael Vazquez <vazquez@im...>  20111103 18:25:39

Hi Xinkang, your mail couldn't be read in the list, I think because of the attached files. I reproduce it below, with the answers to your questions. Regards, Rafa >> Dear GeoPDEs Users, >> I have two questions to ask. >> The first question >> >> I am trying to solve a 3D elasticity problem 'NAFEMS Benchmark >> LE1' by using GeoPDEs. > >> The geometry structure is created by the file 'NAFEMS_LE1.txt' as >> follows: >> # NURBS representation of quarter of a elliptical ring of thickness 0.1, >> # (x/3.25)^2+(y/2.75)^2=1;(x/2)^2+(y/1)^2=1 >> # dimention_of_the_geometry numbers_of_patches >> 3 1 >> # the degree in each Cartesian direction >> 1 2 1 >> # the number of control points in each direction >> 2 3 2 >> # knots >> 0.00000 0.00000 1.00000 1.00000 >> 0.00000 0.00000 0.00000 1.00000 1.00000 1.00000 >> 0.00000 0.00000 1.00000 1.00000 >> # X=x*w >> # y=y*w >> # Z=z*w >> 2.00000 3.25000 1.414213562 2.298097039 0.00000 0.00000 >> 2.00000 3.25000 1.414213562 2.298097039 0.00000 0.00000 >> 0.00000 0.00000 0.707106781 1.944543648 1.00000 2.75000 >> 0.00000 0.00000 0.707106781 1.944543648 1.00000 2.75000 >> 0.00000 0.00000 0.000000000 0.000000000 0.00000 0.00000 0.10000 >> 0.10000 0.0707106781 0.0707106781 0.10000 0.10000 >> # weights >> 1.00000 1.00000 0.707106781 0.707106781 1.00000 1.00000 >> 1.00000 1.00000 0.707106781 0.707106781 1.00000 1.00000 >> Acoording to the structure I built and the model 'NAFEMS Benchmark >> LE1' , we can see that 'problem_data.press_sides=[2]' and >> 'problem_data.symm_sides=[3 4]'; so I use the following code to solve >> this problem: >> % NAFEMS_LE1_3D: solve the linear elasticity problem NAFEMSLE1. >> % 1) PHYSICAL DATA OF THE PROBLEM >> clear problem_data >> % Physical domain, defined as NURBS map given in a text file >> problem_data.geo_name = 'NAFEMS_LE1_3D.txt'; >> % Type of boundary conditions for each side of the domain >> problem_data.nmnn_sides = []; >> problem_data.drchlt_sides = []; >> problem_data.press_sides = [2]; >> problem_data.symm_sides = [3 4]; >> % Physical parameters >> E = 2.1e11; nu = 0.3; >> problem_data.lam = @(x, y, z) ((nu*E)/((1+nu)*(12*nu)) * ones (size >> (x))); >> problem_data.mu = @(x, y, z) (E/(2*(1+nu)) * ones (size (x))); >> % Source and boundary terms >> P = 1.0e7; >> problem_data.f = @(x, y, z) zeros (3, size (x, 1), size (x, 2)); >> problem_data.g = @(x, y, z, ind) zeros (3, size (x, 1), size (x, 2)); >> problem_data.h = @(x, y, z, ind) zeros (3, size (x, 1), size (x, 2)); >> problem_data.p = @(x, y, z, ind) P * ones (size (x)); >> % 2) CHOICE OF THE DISCRETIZATION PARAMETERS >> clear method_data >> method_data.degree = [3 3 3]; % Degree of the bsplines >> method_data.regularity = [2 2 2]; % Regularity of the splines >> method_data.nsub = [10 1 50]; % Number of subdivisions >> method_data.nquad = [4 4 4]; % Points for the Gaussian >> quadrature rule >> % 3) CALL TO THE SOLVER >> [geometry, msh, space, u] = solve_linear_elasticity_3d (problem_data, >> method_data); >> Unfortunately, all the numbers in {u}, which represents the >> displacements of all nodes, are all zero. > >> I can't figure out what's wrong with my code, however, I guess >> someting must be wrong with 'the Source and boundary terms' or >> 'boundary conditions for each side of the domain'. Could you give me >> an advise of how to correct the mistakes? I had not noticed. We implemented the pressure condition for the 2D case, but not for the 3D case. Well, you only have to copy the pressure condition part from "solve_plane_strain_2d" to "solve_linear_elasticity_3d", and adapt it to the 3D case (compute the coordinate z, and make the pressure also a function of z). This should be easy, but if you need any help just ask me. I will add this in the next release. >> The second question >> I have run the example of 'ex_plane_strain_ring' and obtain the >> displacement of all the nodes. How can I find a function to compute >> all the components of stress of all the nodes, and export them to VTK >> format for plotting? The code as it is doesn't have a function to export the components of the stress, but it shouldn't be difficult to implement one by yourself. You will need two functions, similar to sp_eval and sp_eval_msh (or sp_eval_div and sp_eval_div_msh). The main differences should be in the second function (the one with _msh), because instead of using the field shape_functions you will have to use shape_function_gradients to compute the components of the stress. Once you have these two functions, exporting to VTK format should be easy, with functions similar (maybe the same?) to sp_to_vtk and msh_to_vtk. >> Thanks a lot. >> Best regards. >>  >> Xinkang Li >> Zhejiang University >> Hangzhou Zhejiang >> China >> TEL：13456979712 >> Email: lixinkang1987@... <mailto:lixinkang1987@...> >> >> > 