Merged changes in the trunk up to revision 26260.
authorTamito Kajiyama <rd6t-kjym@asahi-net.or.jp>
Mon, 25 Jan 2010 21:06:08 +0000 (21:06 +0000)
committerTamito Kajiyama <rd6t-kjym@asahi-net.or.jp>
Mon, 25 Jan 2010 21:06:08 +0000 (21:06 +0000)
345 files changed:
CMakeLists.txt
config/aix4-config.py [new file with mode: 0644]
config/irix6-config.py
config/linuxcross-config.py
config/openbsd3-config.py
config/sunos5-config.py
intern/audaspace/CMakeLists.txt
intern/audaspace/OpenAL/AUD_OpenALDevice.cpp
intern/ghost/SConscript
intern/ghost/intern/GHOST_WindowX11.cpp
intern/guardedalloc/intern/mallocn.c
intern/smoke/intern/EIGENVALUE_HELPER.cpp [new file with mode: 0644]
intern/smoke/intern/EIGENVALUE_HELPER.h
intern/smoke/intern/FLUID_3D.cpp
intern/smoke/intern/FLUID_3D.h
intern/smoke/intern/FLUID_3D_SOLVERS.cpp
intern/smoke/intern/FLUID_3D_STATIC.cpp
intern/smoke/intern/LU_HELPER.cpp [new file with mode: 0644]
intern/smoke/intern/LU_HELPER.h
intern/smoke/intern/WTURBULENCE.cpp
intern/smoke/intern/smoke_API.cpp
intern/string/intern/STR_String.cpp
projectfiles_vc9/blender/BPY_python/BPY_python.vcproj
projectfiles_vc9/blender/blenlib/BLI_blenlib.vcproj
projectfiles_vc9/blender/editors/ED_editors.vcproj
projectfiles_vc9/blender/nodes/nodes.vcproj
release/datafiles/blenderbuttons
release/scripts/io/engine_render_pov.py
release/scripts/io/export_3ds.py
release/scripts/io/export_fbx.py
release/scripts/io/export_obj.py
release/scripts/io/export_x3d.py
release/scripts/io/import_anim_bvh.py
release/scripts/modules/bpy/__init__.py
release/scripts/modules/bpy/app.py
release/scripts/modules/bpy/ops.py
release/scripts/modules/bpy/utils.py
release/scripts/modules/bpy_types.py
release/scripts/modules/retopo.py
release/scripts/modules/rigify/__init__.py
release/scripts/modules/rigify/arm_biped.py
release/scripts/modules/rigify/delta.py
release/scripts/modules/rigify/eye_balls.py [new file with mode: 0644]
release/scripts/modules/rigify/eye_lid.py [new file with mode: 0644]
release/scripts/modules/rigify/finger_curl.py
release/scripts/modules/rigify/leg_biped.py
release/scripts/modules/rigify/leg_quadruped.py
release/scripts/modules/rigify/mouth.py
release/scripts/modules/rigify/neck_flex.py
release/scripts/modules/rigify/palm_curl.py
release/scripts/modules/rigify/shape_key_distance.py [new file with mode: 0644]
release/scripts/modules/rigify/shape_key_rotdiff.py [new file with mode: 0644]
release/scripts/modules/rigify/spine_pivot_flex.py
release/scripts/modules/rigify_utils.py
release/scripts/modules/rna_info.py
release/scripts/op/mesh.py
release/scripts/op/mesh_skin.py
release/scripts/op/object_randomize_transform.py [new file with mode: 0644]
release/scripts/op/uv.py
release/scripts/op/uvcalc_smart_project.py
release/scripts/templates/gamelogic.py
release/scripts/templates/panel_simple.py [new file with mode: 0644]
release/scripts/ui/properties_animviz.py [new file with mode: 0644]
release/scripts/ui/properties_data_armature.py
release/scripts/ui/properties_data_lamp.py
release/scripts/ui/properties_data_modifier.py
release/scripts/ui/properties_physics_smoke.py
release/scripts/ui/properties_scene.py
release/scripts/ui/properties_texture.py
release/scripts/ui/space_dopesheet.py
release/scripts/ui/space_graph.py
release/scripts/ui/space_image.py
release/scripts/ui/space_info.py
release/scripts/ui/space_nla.py
release/scripts/ui/space_text.py
release/scripts/ui/space_time.py
release/scripts/ui/space_userpref.py
release/scripts/ui/space_view3d.py
release/scripts/ui/space_view3d_toolbar.py
source/Makefile
source/blender/blenkernel/BKE_DerivedMesh.h
source/blender/blenkernel/BKE_anim.h
source/blender/blenkernel/BKE_blender.h
source/blender/blenkernel/BKE_colortools.h
source/blender/blenkernel/BKE_image.h
source/blender/blenkernel/BKE_modifier.h
source/blender/blenkernel/BKE_nla.h
source/blender/blenkernel/BKE_node.h
source/blender/blenkernel/BKE_object.h
source/blender/blenkernel/BKE_sequencer.h
source/blender/blenkernel/BKE_utildefines.h
source/blender/blenkernel/intern/BME_Customdata.c
source/blender/blenkernel/intern/BME_mesh.c
source/blender/blenkernel/intern/DerivedMesh.c
source/blender/blenkernel/intern/anim.c
source/blender/blenkernel/intern/anim_sys.c
source/blender/blenkernel/intern/armature.c
source/blender/blenkernel/intern/blender.c
source/blender/blenkernel/intern/boids.c
source/blender/blenkernel/intern/brush.c
source/blender/blenkernel/intern/colortools.c
source/blender/blenkernel/intern/constraint.c
source/blender/blenkernel/intern/context.c
source/blender/blenkernel/intern/curve.c
source/blender/blenkernel/intern/customdata.c
source/blender/blenkernel/intern/effect.c
source/blender/blenkernel/intern/fcurve.c
source/blender/blenkernel/intern/idprop.c
source/blender/blenkernel/intern/image.c
source/blender/blenkernel/intern/modifier.c
source/blender/blenkernel/intern/multires.c
source/blender/blenkernel/intern/nla.c
source/blender/blenkernel/intern/node.c
source/blender/blenkernel/intern/object.c
source/blender/blenkernel/intern/particle_system.c
source/blender/blenkernel/intern/pointcache.c
source/blender/blenkernel/intern/scene.c
source/blender/blenkernel/intern/smoke.c
source/blender/blenkernel/intern/subsurf_ccg.c
source/blender/blenkernel/intern/unit.c
source/blender/blenlib/BLI_ghash.h
source/blender/blenlib/BLI_gsqueue.h
source/blender/blenlib/BLI_math_base.h
source/blender/blenlib/BLI_math_color.h
source/blender/blenlib/BLI_math_geom.h
source/blender/blenlib/BLI_math_inline.h [new file with mode: 0644]
source/blender/blenlib/BLI_math_vector.h
source/blender/blenlib/BLI_mempool.h
source/blender/blenlib/BLI_threads.h
source/blender/blenlib/intern/BLI_ghash.c
source/blender/blenlib/intern/BLI_kdopbvh.c
source/blender/blenlib/intern/BLI_mempool.c
source/blender/blenlib/intern/gsqueue.c
source/blender/blenlib/intern/math_base.c
source/blender/blenlib/intern/math_base_inline.c [new file with mode: 0644]
source/blender/blenlib/intern/math_color.c
source/blender/blenlib/intern/math_geom.c
source/blender/blenlib/intern/math_vector.c
source/blender/blenlib/intern/math_vector_inline.c
source/blender/blenlib/intern/threads.c
source/blender/blenloader/intern/readfile.c
source/blender/blenloader/intern/writefile.c
source/blender/editors/CMakeLists.txt
source/blender/editors/animation/anim_channels_defines.c
source/blender/editors/animation/anim_channels_edit.c
source/blender/editors/animation/anim_draw.c
source/blender/editors/animation/anim_filter.c
source/blender/editors/animation/anim_ops.c
source/blender/editors/animation/keyframes_draw.c
source/blender/editors/animation/keyframes_edit.c
source/blender/editors/animation/keyframing.c
source/blender/editors/armature/armature_intern.h
source/blender/editors/armature/armature_ops.c
source/blender/editors/armature/editarmature.c
source/blender/editors/armature/poseSlide.c
source/blender/editors/armature/poselib.c
source/blender/editors/armature/poseobject.c
source/blender/editors/curve/editcurve.c
source/blender/editors/datafiles/B.blend.c
source/blender/editors/datafiles/blenderbuttons.c
source/blender/editors/include/ED_anim_api.h
source/blender/editors/include/ED_fileselect.h
source/blender/editors/include/ED_object.h
source/blender/editors/include/UI_icons.h
source/blender/editors/include/UI_interface.h
source/blender/editors/include/UI_resources.h
source/blender/editors/interface/interface.c
source/blender/editors/interface/interface_draw.c
source/blender/editors/interface/interface_handlers.c
source/blender/editors/interface/interface_icons.c
source/blender/editors/interface/interface_intern.h
source/blender/editors/interface/interface_panel.c
source/blender/editors/interface/interface_regions.c
source/blender/editors/interface/interface_style.c
source/blender/editors/interface/interface_templates.c
source/blender/editors/interface/interface_widgets.c
source/blender/editors/interface/resources.c
source/blender/editors/interface/view2d.c
source/blender/editors/mesh/editmesh.c
source/blender/editors/mesh/editmesh_add.c
source/blender/editors/mesh/editmesh_lib.c
source/blender/editors/mesh/editmesh_loop.c
source/blender/editors/object/object_add.c
source/blender/editors/object/object_bake.c
source/blender/editors/object/object_edit.c
source/blender/editors/object/object_intern.h
source/blender/editors/object/object_ops.c
source/blender/editors/object/object_relations.c
source/blender/editors/object/object_vgroup.c
source/blender/editors/physics/particle_edit.c
source/blender/editors/physics/particle_object.c
source/blender/editors/screen/area.c
source/blender/editors/screen/screen_edit.c
source/blender/editors/screen/screen_ops.c
source/blender/editors/screen/screendump.c
source/blender/editors/sculpt_paint/paint_image.c
source/blender/editors/sculpt_paint/paint_vertex.c
source/blender/editors/sculpt_paint/sculpt.c
source/blender/editors/sculpt_paint/sculpt_intern.h
source/blender/editors/space_action/action_edit.c
source/blender/editors/space_action/action_select.c
source/blender/editors/space_action/space_action.c
source/blender/editors/space_api/spacetypes.c
source/blender/editors/space_buttons/space_buttons.c
source/blender/editors/space_file/file_draw.c
source/blender/editors/space_file/file_intern.h
source/blender/editors/space_file/file_ops.c
source/blender/editors/space_file/space_file.c
source/blender/editors/space_graph/graph_buttons.c
source/blender/editors/space_graph/graph_edit.c
source/blender/editors/space_graph/graph_select.c
source/blender/editors/space_graph/space_graph.c
source/blender/editors/space_image/image_buttons.c
source/blender/editors/space_image/image_draw.c
source/blender/editors/space_image/image_intern.h
source/blender/editors/space_image/image_ops.c
source/blender/editors/space_image/space_image.c
source/blender/editors/space_logic/Makefile
source/blender/editors/space_logic/SConscript
source/blender/editors/space_logic/logic_buttons.c
source/blender/editors/space_logic/logic_intern.h
source/blender/editors/space_logic/space_logic.c
source/blender/editors/space_nla/nla_buttons.c
source/blender/editors/space_nla/nla_channels.c
source/blender/editors/space_nla/nla_edit.c
source/blender/editors/space_nla/nla_intern.h
source/blender/editors/space_nla/nla_ops.c
source/blender/editors/space_nla/nla_select.c
source/blender/editors/space_nla/space_nla.c
source/blender/editors/space_node/drawnode.c
source/blender/editors/space_node/node_draw.c
source/blender/editors/space_node/node_edit.c
source/blender/editors/space_node/node_header.c
source/blender/editors/space_node/node_intern.h
source/blender/editors/space_outliner/outliner.c
source/blender/editors/space_outliner/space_outliner.c
source/blender/editors/space_time/space_time.c
source/blender/editors/space_time/time_ops.c
source/blender/editors/space_view3d/drawanimviz.c
source/blender/editors/space_view3d/drawarmature.c
source/blender/editors/space_view3d/drawobject.c
source/blender/editors/space_view3d/space_view3d.c
source/blender/editors/space_view3d/view3d_buttons.c
source/blender/editors/space_view3d/view3d_draw.c
source/blender/editors/space_view3d/view3d_edit.c
source/blender/editors/space_view3d/view3d_intern.h
source/blender/editors/space_view3d/view3d_ops.c
source/blender/editors/space_view3d/view3d_select.c
source/blender/editors/space_view3d/view3d_snap.c
source/blender/editors/space_view3d/view3d_toolbar.c
source/blender/editors/space_view3d/view3d_view.c
source/blender/editors/transform/transform.c
source/blender/editors/transform/transform.h
source/blender/editors/transform/transform_constraints.c
source/blender/editors/transform/transform_conversions.c
source/blender/editors/transform/transform_generics.c
source/blender/editors/transform/transform_numinput.c
source/blender/editors/transform/transform_snap.c
source/blender/ikplugin/intern/iksolver_plugin.c
source/blender/makesdna/DNA_action_types.h
source/blender/makesdna/DNA_anim_types.h
source/blender/makesdna/DNA_brush_types.h
source/blender/makesdna/DNA_color_types.h
source/blender/makesdna/DNA_lamp_types.h
source/blender/makesdna/DNA_node_types.h
source/blender/makesdna/DNA_scene_types.h
source/blender/makesdna/DNA_smoke_types.h
source/blender/makesdna/DNA_space_types.h
source/blender/makesdna/DNA_texture_types.h
source/blender/makesdna/DNA_userdef_types.h
source/blender/makesdna/DNA_view3d_types.h
source/blender/makesdna/intern/SConscript
source/blender/makesrna/RNA_access.h
source/blender/makesrna/RNA_define.h
source/blender/makesrna/RNA_types.h
source/blender/makesrna/intern/CMakeLists.txt
source/blender/makesrna/intern/SConscript
source/blender/makesrna/intern/makesrna.c
source/blender/makesrna/intern/rna_access.c
source/blender/makesrna/intern/rna_action.c
source/blender/makesrna/intern/rna_action_api.c
source/blender/makesrna/intern/rna_brush.c
source/blender/makesrna/intern/rna_camera.c
source/blender/makesrna/intern/rna_color.c
source/blender/makesrna/intern/rna_constraint.c
source/blender/makesrna/intern/rna_curve.c
source/blender/makesrna/intern/rna_define.c
source/blender/makesrna/intern/rna_fcurve.c
source/blender/makesrna/intern/rna_image.c
source/blender/makesrna/intern/rna_internal.h
source/blender/makesrna/intern/rna_lamp.c
source/blender/makesrna/intern/rna_material.c
source/blender/makesrna/intern/rna_modifier.c
source/blender/makesrna/intern/rna_nla.c
source/blender/makesrna/intern/rna_nodetree.c
source/blender/makesrna/intern/rna_nodetree_types.h
source/blender/makesrna/intern/rna_object_api.c
source/blender/makesrna/intern/rna_pose.c
source/blender/makesrna/intern/rna_rna.c
source/blender/makesrna/intern/rna_scene.c
source/blender/makesrna/intern/rna_scene_api.c
source/blender/makesrna/intern/rna_sculpt_paint.c
source/blender/makesrna/intern/rna_smoke.c
source/blender/makesrna/intern/rna_space.c
source/blender/makesrna/intern/rna_texture.c
source/blender/makesrna/intern/rna_timeline.c
source/blender/makesrna/intern/rna_ui_api.c
source/blender/makesrna/intern/rna_userdef.c
source/blender/makesrna/intern/rna_wm.c
source/blender/nodes/CMP_node.h
source/blender/nodes/intern/CMP_nodes/CMP_colorbalance.c [new file with mode: 0644]
source/blender/nodes/intern/CMP_nodes/CMP_huecorrect.c [new file with mode: 0644]
source/blender/nodes/intern/CMP_nodes/CMP_image.c
source/blender/nodes/intern/CMP_nodes/CMP_splitViewer.c
source/blender/nodes/intern/CMP_nodes/CMP_viewer.c
source/blender/nodes/intern/node_util.h
source/blender/python/doc/epy/BGL.py [new file with mode: 0644]
source/blender/python/doc/epy/Geometry.py [new file with mode: 0644]
source/blender/python/doc/epy/IDProp.py [new file with mode: 0644]
source/blender/python/doc/epy/Mathutils.py [new file with mode: 0644]
source/blender/python/doc/epy/testbgl.py [new file with mode: 0644]
source/blender/python/generic/Geometry.c
source/blender/python/generic/Mathutils.c
source/blender/python/generic/euler.c
source/blender/python/generic/matrix.c
source/blender/python/generic/quat.c
source/blender/python/generic/vector.c
source/blender/python/intern/bpy_array.c
source/blender/python/intern/bpy_interface.c
source/blender/python/intern/bpy_operator.c
source/blender/python/intern/bpy_operator_wrap.c
source/blender/python/intern/bpy_props.c [new file with mode: 0644]
source/blender/python/intern/bpy_props.h [new file with mode: 0644]
source/blender/python/intern/bpy_rna.c
source/blender/python/intern/bpy_rna.h
source/blender/python/sphinx_doc_gen.py
source/blender/render/intern/source/texture.c
source/blender/render/intern/source/voxeldata.c
source/blender/windowmanager/WM_types.h
source/blender/windowmanager/intern/wm_event_system.c
source/blender/windowmanager/intern/wm_files.c
source/blender/windowmanager/intern/wm_operators.c
source/blender/windowmanager/intern/wm_window.c
source/gameengine/Converter/BL_BlenderDataConversion.cpp
source/nan_definitions.mk

index 96df995f2e7467541e44bcf3d73d34168534494f..45b13e4b647e114e32710a73a79e7db527050892 100644 (file)
@@ -248,7 +248,7 @@ IF(UNIX AND NOT APPLE)
        SET(PLATFORM_LINKFLAGS "-pthread")
 
        # Better warnings
-       SET(C_WARNINGS "-Wall -Wno-char-subscripts -Wpointer-arith -Wcast-align -Wdeclaration-after-statement")
+       SET(C_WARNINGS "-Wall -Wno-char-subscripts -Wpointer-arith -Wcast-align -Wdeclaration-after-statement -Wno-unknown-pragmas")
        SET(CXX_WARNINGS "-Wall -Wno-invalid-offsetof -Wno-sign-compare")
 
        INCLUDE_DIRECTORIES(${JPEG_INCLUDE_DIR} ${PNG_INCLUDE_DIR} ${ZLIB_INCLUDE_DIR} )
@@ -601,7 +601,7 @@ IF(APPLE)
        ENDIF(CMAKE_OSX_ARCHITECTURES MATCHES "i386")
 
        # Better warnings
-       SET(C_WARNINGS "-Wall -Wno-char-subscripts -Wpointer-arith -Wcast-align -Wdeclaration-after-statement")
+       SET(C_WARNINGS "-Wall -Wno-char-subscripts -Wpointer-arith -Wcast-align -Wdeclaration-after-statement -Wno-unknown-pragmas")
        SET(CXX_WARNINGS "-Wall -Wno-invalid-offsetof -Wno-sign-compare")
 
 ENDIF(APPLE)
diff --git a/config/aix4-config.py b/config/aix4-config.py
new file mode 100644 (file)
index 0000000..3a3db39
--- /dev/null
@@ -0,0 +1,225 @@
+import os
+
+LCGDIR = os.getcwd()+"/../lib/aix-4.3-ppc"
+LIBDIR = LCGDIR
+print LCGDIR
+
+WITH_BF_VERSE = 'false'
+BF_VERSE_INCLUDE = "#extern/verse/dist"
+
+BF_PYTHON = LCGDIR+'/python'
+BF_PYTHON_VERSION = '3.1'
+WITH_BF_STATICPYTHON = 'true'
+BF_PYTHON_INC = '${BF_PYTHON}/include/python${BF_PYTHON_VERSION}'
+BF_PYTHON_BINARY = '${BF_PYTHON}/bin/python${BF_PYTHON_VERSION}'
+BF_PYTHON_LIB = 'python${BF_PYTHON_VERSION}' #BF_PYTHON+'/lib/python'+BF_PYTHON_VERSION+'/config/libpython'+BF_PYTHON_VERSION+'.a'
+BF_PYTHON_LINKFLAGS = ['-Xlinker', '-export-dynamic']
+BF_PYTHON_LIB_STATIC = '${BF_PYTHON}/lib/python2.5/config/libpython${BF_PYTHON_VERSION}.a'
+
+WITH_BF_OPENAL = 'false'
+WITH_BF_STATICOPENAL = 'false'
+BF_OPENAL = LCGDIR+'/openal'
+BF_OPENAL_INC = '${BF_OPENAL}/include'
+BF_OPENAL_LIB = 'openal'
+BF_OPENAL_LIB_STATIC = '${BF_OPENAL}/lib/libopenal.a'
+BF_OPENAL_LIBPATH = LIBDIR + '/lib'
+
+BF_CXX = '/usr'
+WITH_BF_STATICCXX = 'false'
+BF_CXX_LIB_STATIC = '${BF_CXX}/lib/libstdc++.a'
+
+WITH_BF_SDL = 'false'
+BF_SDL = LCGDIR+'/sdl' #$(shell sdl-config --prefix)
+BF_SDL_INC = '${BF_SDL}/include/SDL' #$(shell $(BF_SDL)/bin/sdl-config --cflags)
+BF_SDL_LIB = 'SDL audio iconv charset' #BF_SDL #$(shell $(BF_SDL)/bin/sdl-config --libs) -lSDL_mixer
+BF_SDL_LIBPATH = '${BF_SDL}/lib'
+
+WITH_BF_OPENEXR = 'false'
+WITH_BF_STATICOPENEXR = 'false'
+BF_OPENEXR = '/usr'
+# when compiling with your own openexr lib you might need to set...
+# BF_OPENEXR_INC = '${BF_OPENEXR}/include/OpenEXR ${BF_OPENEXR}/include'
+
+BF_OPENEXR_INC = '${BF_OPENEXR}/include/OpenEXR'
+BF_OPENEXR_LIB = 'Half IlmImf Iex Imath '
+BF_OPENEXR_LIB_STATIC = '${BF_OPENEXR}/lib/libHalf.a ${BF_OPENEXR}/lib/libIlmImf.a ${BF_OPENEXR}/lib/libIex.a ${BF_OPENEXR}/lib/libImath.a ${BF_OPENEXR}/lib/libIlmThread.a'
+# BF_OPENEXR_LIBPATH = '${BF_OPENEXR}/lib'
+
+
+WITH_BF_DDS = 'false'
+
+WITH_BF_JPEG = 'false'
+BF_JPEG = LCGDIR+'/jpeg'
+BF_JPEG_INC = '${BF_JPEG}/include'
+BF_JPEG_LIB = 'jpeg'
+BF_JPEG_LIBPATH = '${BF_JPEG}/lib'
+
+WITH_BF_PNG = 'false'
+BF_PNG = LCGDIR+"/png"
+BF_PNG_INC = '${BF_PNG}/include'
+BF_PNG_LIB = 'png'
+BF_PNG_LIBPATH = '${BF_PNG}/lib'
+
+BF_TIFF = '/usr/nekoware'
+BF_TIFF_INC = '${BF_TIFF}/include'
+
+WITH_BF_ZLIB = 'true'
+BF_ZLIB = LCGDIR+"/zlib"
+BF_ZLIB_INC = '${BF_ZLIB}/include'
+BF_ZLIB_LIB = 'z'
+BF_ZLIB_LIBPATH = '${BF_ZLIB}/lib'
+
+WITH_BF_INTERNATIONAL = 'false'
+
+BF_GETTEXT = LCGDIR+'/gettext'
+BF_GETTEXT_INC = '${BF_GETTEXT}/include'
+BF_GETTEXT_LIB = 'gettextpo intl'
+BF_GETTEXT_LIBPATH = '${BF_GETTEXT}/lib'
+
+WITH_BF_FTGL = 'false'
+BF_FTGL = '#extern/bFTGL'
+BF_FTGL_INC = '${BF_FTGL}/include'
+BF_FTGL_LIB = 'extern_ftgl'
+
+WITH_BF_GAMEENGINE='false'
+
+WITH_BF_ODE = 'false'
+BF_ODE = LIBDIR + '/ode'
+BF_ODE_INC = BF_ODE + '/include'
+BF_ODE_LIB = BF_ODE + '/lib/libode.a'
+
+WITH_BF_BULLET = 'true'
+BF_BULLET = '#extern/bullet2/src'
+BF_BULLET_INC = '${BF_BULLET}'
+BF_BULLET_LIB = 'extern_bullet'
+
+BF_SOLID = '#extern/solid'
+BF_SOLID_INC = '${BF_SOLID}'
+BF_SOLID_LIB = 'extern_solid'
+
+WITH_BF_YAFRAY = 'true'
+
+#WITH_BF_NSPR = 'true'
+#BF_NSPR = $(LIBDIR)/nspr
+#BF_NSPR_INC = -I$(BF_NSPR)/include -I$(BF_NSPR)/include/nspr
+#BF_NSPR_LIB = 
+
+# Uncomment the following line to use Mozilla inplace of netscape
+#CPPFLAGS += -DMOZ_NOT_NET
+# Location of MOZILLA/Netscape header files...
+#BF_MOZILLA = $(LIBDIR)/mozilla
+#BF_MOZILLA_INC = -I$(BF_MOZILLA)/include/mozilla/nspr -I$(BF_MOZILLA)/include/mozilla -I$(BF_MOZILLA)/include/mozilla/xpcom -I$(BF_MOZILLA)/include/mozilla/idl
+#BF_MOZILLA_LIB =
+# Will fall back to look in BF_MOZILLA_INC/nspr and BF_MOZILLA_LIB
+# if this is not set.
+#
+# Be paranoid regarding library creation (do not update archives)
+#BF_PARANOID = 'true'
+
+# enable freetype2 support for text objects
+BF_FREETYPE = LCGDIR+'/freetype'
+BF_FREETYPE_INC = '${BF_FREETYPE}/include ${BF_FREETYPE}/include/freetype2'
+BF_FREETYPE_LIB = 'freetype'
+BF_FREETYPE_LIBPATH = '${BF_FREETYPE}/lib'
+
+WITH_BF_QUICKTIME = 'false' # -DWITH_QUICKTIME
+BF_QUICKTIME = '/usr/local'
+BF_QUICKTIME_INC = '${BF_QUICKTIME}/include'
+
+WITH_BF_ICONV = 'false'
+BF_ICONV = LIBDIR + "/iconv"
+BF_ICONV_INC = '${BF_ICONV}/include'
+BF_ICONV_LIB = 'iconv charset'
+BF_ICONV_LIBPATH = '${BF_ICONV}/lib'
+
+WITH_BF_BINRELOC = 'false'
+
+# enable ffmpeg  support
+WITH_BF_FFMPEG = 'false'  # -DWITH_FFMPEG
+# Uncomment the following two lines to use system's ffmpeg
+BF_FFMPEG = LCGDIR+'/ffmpeg'
+BF_FFMPEG_LIB = 'avformat avcodec swscale avutil avdevice faad faac vorbis x264 ogg mp3lame z'
+BF_FFMPEG_INC = '${BF_FFMPEG}/include'
+BF_FFMPEG_LIBPATH='${BF_FFMPEG}/lib'
+
+# enable ogg, vorbis and theora in ffmpeg
+WITH_BF_OGG = 'false'  # -DWITH_OGG 
+BF_OGG = '/usr'
+BF_OGG_INC = '${BF_OGG}/include'
+BF_OGG_LIB = 'ogg vorbis theoraenc theoradec'
+
+WITH_BF_OPENJPEG = 'false' 
+BF_OPENJPEG = '#extern/libopenjpeg'
+BF_OPENJPEG_LIB = ''
+BF_OPENJPEG_INC = '${BF_OPENJPEG}'
+BF_OPENJPEG_LIBPATH='${BF_OPENJPEG}/lib'
+
+WITH_BF_REDCODE = 'false'  
+BF_REDCODE = '#extern/libredcode'
+BF_REDCODE_LIB = ''
+BF_REDCODE_INC = '${BF_REDCODE}/include'
+BF_REDCODE_LIBPATH='${BF_REDCODE}/lib'
+
+# Mesa Libs should go here if your using them as well....
+WITH_BF_STATICOPENGL = 'false'
+BF_OPENGL = '/usr'
+BF_OPENGL_INC = '${BF_OPENGL}/include'
+BF_OPENGL_LIB = 'GL GLU X11 Xi Xext'
+BF_OPENGL_LIBPATH = '/usr/X11R6/lib'
+BF_OPENGL_LIB_STATIC = '${BF_OPENGL}/libGL.a ${BF_OPENGL}/libGLU.a ${BF_OPENGL}/libXxf86vm.a ${BF_OPENGL}/libX11.a ${BF_OPENGL}/libXi.a ${BF_OPENGL}/libXext.a ${BF_OPENGL}/libXxf86vm.a'
+
+
+CC = 'gcc'
+CXX = 'g++'
+
+CCFLAGS = [ '-pipe', '-funsigned-char', '-fno-strict-aliasing' ]
+
+CPPFLAGS = [ '-DXP_UNIX', '-DWIN32', '-DFREE_WINDOWS' ]
+CXXFLAGS = ['-pipe', '-funsigned-char', '-fno-strict-aliasing' ]
+REL_CFLAGS = [ '-O2' ]
+REL_CCFLAGS = [ '-O2' ]
+C_WARN = [ '-Wall' , '-Wno-char-subscripts', '-Wdeclaration-after-statement' ]
+
+CC_WARN = [ '-Wall' ]
+
+
+
+##BF_DEPEND = 'true'
+##
+##AR = ar
+##ARFLAGS = ruv
+##ARFLAGSQUIET = ru
+##
+
+##FIX_STUBS_WARNINGS = -Wno-unused
+
+LLIBS = 'c m dl pthread dmedia movie'
+##LOPTS = --dynamic
+##DYNLDFLAGS = -shared $(LDFLAGS)
+
+BF_PROFILE_FLAGS = ['-pg','-g']
+BF_PROFILE = 'false'
+
+BF_DEBUG = 'false'
+BF_DEBUG_FLAGS = '-g'
+
+BF_BUILDDIR = '../build/aix4'
+BF_INSTALLDIR='../install/aix4'
+BF_DOCDIR='../install/doc'
+
+#Link against pthread
+LDIRS = []
+LDIRS.append(BF_FREETYPE_LIBPATH)
+LDIRS.append(BF_PNG_LIBPATH)
+LDIRS.append(BF_ZLIB_LIBPATH)
+LDIRS.append(BF_SDL_LIBPATH)
+LDIRS.append(BF_OPENAL_LIBPATH)
+LDIRS.append(BF_ICONV_LIBPATH)
+
+PLATFORM_LINKFLAGS = []
+for x in LDIRS:
+    PLATFORM_LINKFLAGS.append("-L"+x)
+    
+PLATFORM_LINKFLAGS += ['-L${LCGDIR}/jpeg/lib' , '-L/usr/lib32',  '-n32', '-v', '-no_prelink']
+print PLATFORM_LINKFLAGS
+LINKFLAGS= PLATFORM_LINKFLAGS
index ab2ef02c977c65f0be1ef1cfe576b5e6dc653294..df18e0b511f6584208f78cb0310ed40c75d597e4 100644 (file)
@@ -4,7 +4,7 @@ LCGDIR = os.getcwd()+"/../lib/irix-6.5-mips"
 LIBDIR = LCGDIR
 
 BF_PYTHON = LCGDIR+'/python'
-BF_PYTHON_VERSION = '2.5'
+BF_PYTHON_VERSION = '3.1'
 WITH_BF_STATICPYTHON = 'true'
 BF_PYTHON_INC = '${BF_PYTHON}/include/python${BF_PYTHON_VERSION}'
 BF_PYTHON_BINARY = '${BF_PYTHON}/bin/python${BF_PYTHON_VERSION}'
index 3cfa11587004915ae8e22c0464dc8b32a8cc9a6b..571d644a9c4a67824ae91f9c4ca3c26e77ccdece 100644 (file)
@@ -2,7 +2,7 @@ LCGDIR = '#../lib/windows'
 LIBDIR = '${LCGDIR}'
 
 BF_PYTHON = LIBDIR + '/python'
-BF_PYTHON_VERSION = '2.5'
+BF_PYTHON_VERSION = '3.1'
 BF_PYTHON_INC = '${BF_PYTHON}/include/python${BF_PYTHON_VERSION}'
 BF_PYTHON_BINARY = 'python'
 BF_PYTHON_LIB = 'python25'
index 353d30f50b3e90d41310ce33b1f8f4a53b604ad4..0ef9ba5d0a435d2b8b1aa7f12e8e9a7ace1fc6d5 100644 (file)
@@ -2,7 +2,7 @@ LCGDIR = '../lib/openbsd3'
 LIBDIR = '${LCGDIR}'
 
 BF_PYTHON = '/usr/local'
-BF_PYTHON_VERSION = '2.5'
+BF_PYTHON_VERSION = '3.1'
 BF_PYTHON_INC = '${BF_PYTHON}/include/python${BF_PYTHON_VERSION}'
 BF_PYTHON_BINARY = '${BF_PYTHON}/bin/python${BF_PYTHON_VERSION}'
 BF_PYTHON_LIB = 'python${BF_PYTHON_VERSION}'
index 88dce927db4f9240b1a40f8ca4b54774eb1af341..a0713735a5b5fab67ea18657c2335d7658371064 100644 (file)
@@ -2,7 +2,7 @@ LCGDIR = '../lib/sunos5'
 LIBDIR = '${LCGDIR}'
 
 BF_PYTHON = '/usr/local'
-BF_PYTHON_VERSION = '2.5'
+BF_PYTHON_VERSION = '3.1'
 BF_PYTHON_INC = '${BF_PYTHON}/include/python${BF_PYTHON_VERSION}'
 BF_PYTHON_BINARY = '${BF_PYTHON}/bin/python${BF_PYTHON_VERSION}'
 BF_PYTHON_LIB = 'python${BF_PYTHON_VERSION}' #BF_PYTHON+'/lib/python'+BF_PYTHON_VERSION+'/config/libpython'+BF_PYTHON_VERSION+'.a'
index 7a3cb88db9cd42388d9554d31c7e2deb21a51228..0965a467201172cbfe5930306c29a68d4ee494c8 100644 (file)
@@ -20,7 +20,7 @@
 #
 # ***** END LGPL LICENSE BLOCK *****
 
-SET(INC . intern FX SRC ${PTHREADS_INC} ${LIBSAMPLERATE_INC} ${OPENAL_INCLUDE_DIR})
+SET(INC . intern FX SRC ${PTHREADS_INC} ${LIBSAMPLERATE_INC})
 
 FILE(GLOB SRC intern/*.cpp intern/*.h FX/*.cpp SRC/*.cpp)
 
index 087c7849a2c44aa27dba7dcecc11e01a3ebd04c0..0b3e86eda56fb7c4ad91bd0b9eec383626cd3b3c 100644 (file)
@@ -284,9 +284,7 @@ AUD_OpenALDevice::AUD_OpenALDevice(AUD_DeviceSpecs specs, int buffersize)
        }
 #endif
 
-       m_device = alcOpenDevice("ALSA Software");
-       if(m_device == NULL)
-               m_device = alcOpenDevice(NULL);
+       m_device = alcOpenDevice(NULL);
 
        if(!m_device)
                AUD_THROW(AUD_ERROR_OPENAL);
index ff98981a35e6ae2199bf1f276ae7923f2a32aeeb..6713ded0afa85b295596e44c3041bcfd98a40f78 100644 (file)
@@ -14,7 +14,7 @@ if window_system == 'darwin':
 pf = ['GHOST_DisplayManager', 'GHOST_System', 'GHOST_Window', 'GHOST_DropTarget']
 defs=['_USE_MATH_DEFINES']
 
-if window_system in ('linux2', 'openbsd3', 'sunos5', 'freebsd6', 'irix6'):
+if window_system in ('linux2', 'openbsd3', 'sunos5', 'freebsd6', 'irix6', 'aix4', 'aix5'):
        for f in pf:
                try:
                        sources.remove('intern' + os.sep + f + 'Win32.cpp')
index efde9005fe4eba973b74221c667c27c22f97e5f2..6c15d5d4e32a1f33a21cd0424ae13fc7ce84cb47 100644 (file)
@@ -35,7 +35,7 @@
 #include <X11/cursorfont.h>
 #include <X11/Xatom.h>
 
-#if defined(__sun__) || defined( __sun ) || defined (__sparc) || defined (__sparc__)
+#if defined(__sun__) || defined( __sun ) || defined (__sparc) || defined (__sparc__) || defined (_AIX)
 #include <strings.h>
 #endif
 
@@ -1435,7 +1435,7 @@ setWindowCursorGrab(
                                setWindowCursorVisibility(false);
 
                }
-               XGrabPointer(m_display, m_window, True, ButtonPressMask| ButtonReleaseMask|PointerMotionMask, GrabModeAsync, GrabModeAsync, None, None, CurrentTime);
+               XGrabPointer(m_display, m_window, False, ButtonPressMask| ButtonReleaseMask|PointerMotionMask, GrabModeAsync, GrabModeAsync, None, None, CurrentTime);
        }
        else {
                if (m_cursorGrab==GHOST_kGrabHide) {
index d1114af2437c09e7a5a27e41cdd7ade2574701de..446daea1c0389639bb5923c66aefd04a7eb896e4 100644 (file)
@@ -487,8 +487,10 @@ short MEM_testN(void *vmemh) {
        if (membl) membl = MEMNEXT(membl);
 
        while(membl) {
-               if (vmemh == membl+1)
+               if (vmemh == membl+1) {
+                       mem_unlock_thread();
                        return 1;
+               }
 
                if(membl->next)
                        membl= MEMNEXT(membl->next);
@@ -628,8 +630,8 @@ static void rem_memblock(MemHead *memh)
         mmap_in_use -= memh->len;
         if (munmap(memh, memh->len + sizeof(MemHead) + sizeof(MemTail)))
             printf("Couldn't unmap memory %s\n", memh->name);
-    }
-    else {
+    }  
+       else {
                if(malloc_debug_memset && memh->len)
                        memset(memh+1, 255, memh->len);
         free(memh);
diff --git a/intern/smoke/intern/EIGENVALUE_HELPER.cpp b/intern/smoke/intern/EIGENVALUE_HELPER.cpp
new file mode 100644 (file)
index 0000000..f11b579
--- /dev/null
@@ -0,0 +1,885 @@
+
+#include "EIGENVALUE_HELPER.h"
+
+
+void Eigentred2(sEigenvalue& eval) {
+
+   //  This is derived from the Algol procedures tred2 by
+   //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+   //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+   //  Fortran subroutine in EISPACK.
+
+         int n=eval.n;
+
+      for (int j = 0; j < n; j++) {
+         eval.d[j] = eval.V[n-1][j];
+      }
+
+      // Householder reduction to tridiagonal form.
+   
+      for (int i = n-1; i > 0; i--) {
+   
+         // Scale to avoid under/overflow.
+   
+         float scale = 0.0;
+         float h = 0.0;
+         for (int k = 0; k < i; k++) {
+            scale = scale + fabs(eval.d[k]);
+         }
+         if (scale == 0.0) {
+            eval.e[i] = eval.d[i-1];
+            for (int j = 0; j < i; j++) {
+               eval.d[j] = eval.V[i-1][j];
+               eval.V[i][j] = 0.0;
+               eval.V[j][i] = 0.0;
+            }
+         } else {
+   
+            // Generate Householder vector.
+   
+            for (int k = 0; k < i; k++) {
+               eval.d[k] /= scale;
+               h += eval.d[k] * eval.d[k];
+            }
+            float f = eval.d[i-1];
+            float g = sqrt(h);
+            if (f > 0) {
+               g = -g;
+            }
+            eval.e[i] = scale * g;
+            h = h - f * g;
+            eval.d[i-1] = f - g;
+            for (int j = 0; j < i; j++) {
+               eval.e[j] = 0.0;
+            }
+   
+            // Apply similarity transformation to remaining columns.
+   
+            for (int j = 0; j < i; j++) {
+               f = eval.d[j];
+               eval.V[j][i] = f;
+               g = eval.e[j] + eval.V[j][j] * f;
+               for (int k = j+1; k <= i-1; k++) {
+                  g += eval.V[k][j] * eval.d[k];
+                  eval.e[k] += eval.V[k][j] * f;
+               }
+               eval.e[j] = g;
+            }
+            f = 0.0;
+            for (int j = 0; j < i; j++) {
+               eval.e[j] /= h;
+               f += eval.e[j] * eval.d[j];
+            }
+            float hh = f / (h + h);
+            for (int j = 0; j < i; j++) {
+               eval.e[j] -= hh * eval.d[j];
+            }
+            for (int j = 0; j < i; j++) {
+               f = eval.d[j];
+               g = eval.e[j];
+               for (int k = j; k <= i-1; k++) {
+                  eval.V[k][j] -= (f * eval.e[k] + g * eval.d[k]);
+               }
+               eval.d[j] = eval.V[i-1][j];
+               eval.V[i][j] = 0.0;
+            }
+         }
+         eval.d[i] = h;
+      }
+   
+      // Accumulate transformations.
+   
+      for (int i = 0; i < n-1; i++) {
+         eval.V[n-1][i] = eval.V[i][i];
+         eval.V[i][i] = 1.0;
+         float h = eval.d[i+1];
+         if (h != 0.0) {
+            for (int k = 0; k <= i; k++) {
+               eval.d[k] = eval.V[k][i+1] / h;
+            }
+            for (int j = 0; j <= i; j++) {
+               float g = 0.0;
+               for (int k = 0; k <= i; k++) {
+                  g += eval.V[k][i+1] * eval.V[k][j];
+               }
+               for (int k = 0; k <= i; k++) {
+                  eval.V[k][j] -= g * eval.d[k];
+               }
+            }
+         }
+         for (int k = 0; k <= i; k++) {
+            eval.V[k][i+1] = 0.0;
+         }
+      }
+      for (int j = 0; j < n; j++) {
+         eval.d[j] = eval.V[n-1][j];
+         eval.V[n-1][j] = 0.0;
+      }
+      eval.V[n-1][n-1] = 1.0;
+      eval.e[0] = 0.0;
+}
+
+void Eigencdiv(sEigenvalue& eval, float xr, float xi, float yr, float yi) {
+      float r,d;
+      if (fabs(yr) > fabs(yi)) {
+         r = yi/yr;
+         d = yr + r*yi;
+         eval.cdivr = (xr + r*xi)/d;
+         eval.cdivi = (xi - r*xr)/d;
+      } else {
+         r = yr/yi;
+         d = yi + r*yr;
+         eval.cdivr = (r*xr + xi)/d;
+         eval.cdivi = (r*xi - xr)/d;
+      }
+   }
+
+void Eigentql2 (sEigenvalue& eval) {
+
+   //  This is derived from the Algol procedures tql2, by
+   //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+   //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+   //  Fortran subroutine in EISPACK.
+   
+         int n=eval.n;
+
+      for (int i = 1; i < n; i++) {
+         eval.e[i-1] = eval.e[i];
+      }
+      eval.e[n-1] = 0.0;
+   
+      float f = 0.0;
+      float tst1 = 0.0;
+      float eps = pow(2.0,-52.0);
+      for (int l = 0; l < n; l++) {
+
+         // Find small subdiagonal element
+   
+         tst1 = max(tst1,fabs(eval.d[l]) + fabs(eval.e[l]));
+         int m = l;
+
+        // Original while-loop from Java code
+         while (m < n) {
+            if (fabs(eval.e[m]) <= eps*tst1) {
+               break;
+            }
+            m++;
+         }
+
+   
+         // If m == l, d[l] is an eigenvalue,
+         // otherwise, iterate.
+   
+         if (m > l) {
+            int iter = 0;
+            do {
+               iter = iter + 1;  // (Could check iteration count here.)
+   
+               // Compute implicit shift
+
+               float g = eval.d[l];
+               float p = (eval.d[l+1] - g) / (2.0 * eval.e[l]);
+               float r = hypot(p,1.0);
+               if (p < 0) {
+                  r = -r;
+               }
+               eval.d[l] = eval.e[l] / (p + r);
+               eval.d[l+1] = eval.e[l] * (p + r);
+               float dl1 = eval.d[l+1];
+               float h = g - eval.d[l];
+               for (int i = l+2; i < n; i++) {
+                  eval.d[i] -= h;
+               }
+               f = f + h;
+   
+               // Implicit QL transformation.
+   
+               p = eval.d[m];
+               float c = 1.0;
+               float c2 = c;
+               float c3 = c;
+               float el1 = eval.e[l+1];
+               float s = 0.0;
+               float s2 = 0.0;
+               for (int i = m-1; i >= l; i--) {
+                  c3 = c2;
+                  c2 = c;
+                  s2 = s;
+                  g = c * eval.e[i];
+                  h = c * p;
+                  r = hypot(p,eval.e[i]);
+                  eval.e[i+1] = s * r;
+                  s = eval.e[i] / r;
+                  c = p / r;
+                  p = c * eval.d[i] - s * g;
+                  eval.d[i+1] = h + s * (c * g + s * eval.d[i]);
+   
+                  // Accumulate transformation.
+   
+                  for (int k = 0; k < n; k++) {
+                     h = eval.V[k][i+1];
+                     eval.V[k][i+1] = s * eval.V[k][i] + c * h;
+                     eval.V[k][i] = c * eval.V[k][i] - s * h;
+                  }
+               }
+               p = -s * s2 * c3 * el1 * eval.e[l] / dl1;
+               eval.e[l] = s * p;
+               eval.d[l] = c * p;
+   
+               // Check for convergence.
+   
+            } while (fabs(eval.e[l]) > eps*tst1);
+         }
+         eval.d[l] = eval.d[l] + f;
+         eval.e[l] = 0.0;
+      }
+     
+      // Sort eigenvalues and corresponding vectors.
+   
+      for (int i = 0; i < n-1; i++) {
+         int k = i;
+         float p = eval.d[i];
+         for (int j = i+1; j < n; j++) {
+            if (eval.d[j] < p) {
+               k = j;
+               p = eval.d[j];
+            }
+         }
+         if (k != i) {
+            eval.d[k] = eval.d[i];
+            eval.d[i] = p;
+            for (int j = 0; j < n; j++) {
+               p = eval.V[j][i];
+               eval.V[j][i] = eval.V[j][k];
+               eval.V[j][k] = p;
+            }
+         }
+      }
+}
+
+void Eigenorthes (sEigenvalue& eval) {
+   
+      //  This is derived from the Algol procedures orthes and ortran,
+      //  by Martin and Wilkinson, Handbook for Auto. Comp.,
+      //  Vol.ii-Linear Algebra, and the corresponding
+      //  Fortran subroutines in EISPACK.
+   
+         int n=eval.n;
+
+      int low = 0;
+      int high = n-1;
+   
+      for (int m = low+1; m <= high-1; m++) {
+   
+         // Scale column.
+   
+         float scale = 0.0;
+         for (int i = m; i <= high; i++) {
+            scale = scale + fabs(eval.H[i][m-1]);
+         }
+         if (scale != 0.0) {
+   
+            // Compute Householder transformation.
+   
+            float h = 0.0;
+            for (int i = high; i >= m; i--) {
+               eval.ort[i] = eval.H[i][m-1]/scale;
+               h += eval.ort[i] * eval.ort[i];
+            }
+            float g = sqrt(h);
+            if (eval.ort[m] > 0) {
+               g = -g;
+            }
+            h = h - eval.ort[m] * g;
+            eval.ort[m] = eval.ort[m] - g;
+   
+            // Apply Householder similarity transformation
+            // H = (I-u*u'/h)*H*(I-u*u')/h)
+   
+            for (int j = m; j < n; j++) {
+               float f = 0.0;
+               for (int i = high; i >= m; i--) {
+                  f += eval.ort[i]*eval.H[i][j];
+               }
+               f = f/h;
+               for (int i = m; i <= high; i++) {
+                  eval.H[i][j] -= f*eval.ort[i];
+               }
+           }
+   
+           for (int i = 0; i <= high; i++) {
+               float f = 0.0;
+               for (int j = high; j >= m; j--) {
+                  f += eval.ort[j]*eval.H[i][j];
+               }
+               f = f/h;
+               for (int j = m; j <= high; j++) {
+                  eval.H[i][j] -= f*eval.ort[j];
+               }
+            }
+            eval.ort[m] = scale*eval.ort[m];
+            eval.H[m][m-1] = scale*g;
+         }
+      }
+   
+      // Accumulate transformations (Algol's ortran).
+
+      for (int i = 0; i < n; i++) {
+         for (int j = 0; j < n; j++) {
+            eval.V[i][j] = (i == j ? 1.0 : 0.0);
+         }
+      }
+
+      for (int m = high-1; m >= low+1; m--) {
+         if (eval.H[m][m-1] != 0.0) {
+            for (int i = m+1; i <= high; i++) {
+               eval.ort[i] = eval.H[i][m-1];
+            }
+            for (int j = m; j <= high; j++) {
+               float g = 0.0;
+               for (int i = m; i <= high; i++) {
+                  g += eval.ort[i] * eval.V[i][j];
+               }
+               // Double division avoids possible underflow
+               g = (g / eval.ort[m]) / eval.H[m][m-1];
+               for (int i = m; i <= high; i++) {
+                  eval.V[i][j] += g * eval.ort[i];
+               }
+            }
+         }
+      }
+   }
+
+void Eigenhqr2 (sEigenvalue& eval) {
+   
+      //  This is derived from the Algol procedure hqr2,
+      //  by Martin and Wilkinson, Handbook for Auto. Comp.,
+      //  Vol.ii-Linear Algebra, and the corresponding
+      //  Fortran subroutine in EISPACK.
+   
+      // Initialize
+   
+      int nn = eval.n;
+      int n = nn-1;
+      int low = 0;
+      int high = nn-1;
+      float eps = pow(2.0,-52.0);
+      float exshift = 0.0;
+      float p=0,q=0,r=0,s=0,z=0,t,w,x,y;
+   
+      // Store roots isolated by balanc and compute matrix norm
+   
+      float norm = 0.0;
+      for (int i = 0; i < nn; i++) {
+         if ((i < low) || (i > high)) {
+            eval.d[i] = eval.H[i][i];
+            eval.e[i] = 0.0;
+         }
+         for (int j = max(i-1,0); j < nn; j++) {
+            norm = norm + fabs(eval.H[i][j]);
+         }
+      }
+   
+      // Outer loop over eigenvalue index
+   
+      int iter = 0;
+               int totIter = 0;
+      while (n >= low) {
+
+                       // NT limit no. of iterations
+                       totIter++;
+                       if(totIter>100) {
+                               //if(totIter>15) std::cout<<"!!!!iter ABORT !!!!!!! "<<totIter<<"\n"; 
+                               // NT hack/fix, return large eigenvalues
+                               for (int i = 0; i < nn; i++) {
+                                       eval.d[i] = 10000.;
+                                       eval.e[i] = 10000.;
+                               }
+                               return;
+                       }
+   
+         // Look for single small sub-diagonal element
+   
+         int l = n;
+         while (l > low) {
+            s = fabs(eval.H[l-1][l-1]) + fabs(eval.H[l][l]);
+            if (s == 0.0) {
+               s = norm;
+            }
+            if (fabs(eval.H[l][l-1]) < eps * s) {
+               break;
+            }
+            l--;
+         }
+       
+         // Check for convergence
+         // One root found
+   
+         if (l == n) {
+            eval.H[n][n] = eval.H[n][n] + exshift;
+            eval.d[n] = eval.H[n][n];
+            eval.e[n] = 0.0;
+            n--;
+            iter = 0;
+   
+         // Two roots found
+   
+         } else if (l == n-1) {
+            w = eval.H[n][n-1] * eval.H[n-1][n];
+            p = (eval.H[n-1][n-1] - eval.H[n][n]) / 2.0;
+            q = p * p + w;
+            z = sqrt(fabs(q));
+            eval.H[n][n] = eval.H[n][n] + exshift;
+            eval.H[n-1][n-1] = eval.H[n-1][n-1] + exshift;
+            x = eval.H[n][n];
+   
+            // float pair
+   
+            if (q >= 0) {
+               if (p >= 0) {
+                  z = p + z;
+               } else {
+                  z = p - z;
+               }
+               eval.d[n-1] = x + z;
+               eval.d[n] = eval.d[n-1];
+               if (z != 0.0) {
+                  eval.d[n] = x - w / z;
+               }
+               eval.e[n-1] = 0.0;
+               eval.e[n] = 0.0;
+               x = eval.H[n][n-1];
+               s = fabs(x) + fabs(z);
+               p = x / s;
+               q = z / s;
+               r = sqrt(p * p+q * q);
+               p = p / r;
+               q = q / r;
+   
+               // Row modification
+   
+               for (int j = n-1; j < nn; j++) {
+                  z = eval.H[n-1][j];
+                  eval.H[n-1][j] = q * z + p * eval.H[n][j];
+                  eval.H[n][j] = q * eval.H[n][j] - p * z;
+               }
+   
+               // Column modification
+   
+               for (int i = 0; i <= n; i++) {
+                  z = eval.H[i][n-1];
+                  eval.H[i][n-1] = q * z + p * eval.H[i][n];
+                  eval.H[i][n] = q * eval.H[i][n] - p * z;
+               }
+   
+               // Accumulate transformations
+   
+               for (int i = low; i <= high; i++) {
+                  z = eval.V[i][n-1];
+                  eval.V[i][n-1] = q * z + p * eval.V[i][n];
+                  eval.V[i][n] = q * eval.V[i][n] - p * z;
+               }
+   
+            // Complex pair
+   
+            } else {
+               eval.d[n-1] = x + p;
+               eval.d[n] = x + p;
+               eval.e[n-1] = z;
+               eval.e[n] = -z;
+            }
+            n = n - 2;
+            iter = 0;
+   
+         // No convergence yet
+   
+         } else {
+   
+            // Form shift
+   
+            x = eval.H[n][n];
+            y = 0.0;
+            w = 0.0;
+            if (l < n) {
+               y = eval.H[n-1][n-1];
+               w = eval.H[n][n-1] * eval.H[n-1][n];
+            }
+   
+            // Wilkinson's original ad hoc shift
+   
+            if (iter == 10) {
+               exshift += x;
+               for (int i = low; i <= n; i++) {
+                  eval.H[i][i] -= x;
+               }
+               s = fabs(eval.H[n][n-1]) + fabs(eval.H[n-1][n-2]);
+               x = y = 0.75 * s;
+               w = -0.4375 * s * s;
+            }
+
+            // MATLAB's new ad hoc shift
+
+            if (iter == 30) {
+                s = (y - x) / 2.0;
+                s = s * s + w;
+                if (s > 0) {
+                    s = sqrt(s);
+                    if (y < x) {
+                       s = -s;
+                    }
+                    s = x - w / ((y - x) / 2.0 + s);
+                    for (int i = low; i <= n; i++) {
+                       eval.H[i][i] -= s;
+                    }
+                    exshift += s;
+                    x = y = w = 0.964;
+                }
+            }
+   
+            iter = iter + 1;   // (Could check iteration count here.)
+   
+            // Look for two consecutive small sub-diagonal elements
+   
+            int m = n-2;
+            while (m >= l) {
+               z = eval.H[m][m];
+               r = x - z;
+               s = y - z;
+               p = (r * s - w) / eval.H[m+1][m] + eval.H[m][m+1];
+               q = eval.H[m+1][m+1] - z - r - s;
+               r = eval.H[m+2][m+1];
+               s = fabs(p) + fabs(q) + fabs(r);
+               p = p / s;
+               q = q / s;
+               r = r / s;
+               if (m == l) {
+                  break;
+               }
+               if (fabs(eval.H[m][m-1]) * (fabs(q) + fabs(r)) <
+                  eps * (fabs(p) * (fabs(eval.H[m-1][m-1]) + fabs(z) +
+                  fabs(eval.H[m+1][m+1])))) {
+                     break;
+               }
+               m--;
+            }
+   
+            for (int i = m+2; i <= n; i++) {
+               eval.H[i][i-2] = 0.0;
+               if (i > m+2) {
+                  eval.H[i][i-3] = 0.0;
+               }
+            }
+   
+            // Double QR step involving rows l:n and columns m:n
+   
+            for (int k = m; k <= n-1; k++) {
+               int notlast = (k != n-1);
+               if (k != m) {
+                  p = eval.H[k][k-1];
+                  q = eval.H[k+1][k-1];
+                  r = (notlast ? eval.H[k+2][k-1] : 0.0);
+                  x = fabs(p) + fabs(q) + fabs(r);
+                  if (x != 0.0) {
+                     p = p / x;
+                     q = q / x;
+                     r = r / x;
+                  }
+               }
+               if (x == 0.0) {
+                  break;
+               }
+               s = sqrt(p * p + q * q + r * r);
+               if (p < 0) {
+                  s = -s;
+               }
+               if (s != 0) {
+                  if (k != m) {
+                     eval.H[k][k-1] = -s * x;
+                  } else if (l != m) {
+                     eval.H[k][k-1] = -eval.H[k][k-1];
+                  }
+                  p = p + s;
+                  x = p / s;
+                  y = q / s;
+                  z = r / s;
+                  q = q / p;
+                  r = r / p;
+   
+                  // Row modification
+   
+                  for (int j = k; j < nn; j++) {
+                     p = eval.H[k][j] + q * eval.H[k+1][j];
+                     if (notlast) {
+                        p = p + r * eval.H[k+2][j];
+                        eval.H[k+2][j] = eval.H[k+2][j] - p * z;
+                     }
+                     eval.H[k][j] = eval.H[k][j] - p * x;
+                     eval.H[k+1][j] = eval.H[k+1][j] - p * y;
+                  }
+   
+                  // Column modification
+   
+                  for (int i = 0; i <= min(n,k+3); i++) {
+                     p = x * eval.H[i][k] + y * eval.H[i][k+1];
+                     if (notlast) {
+                        p = p + z * eval.H[i][k+2];
+                        eval.H[i][k+2] = eval.H[i][k+2] - p * r;
+                     }
+                     eval.H[i][k] = eval.H[i][k] - p;
+                     eval.H[i][k+1] = eval.H[i][k+1] - p * q;
+                  }
+   
+                  // Accumulate transformations
+   
+                  for (int i = low; i <= high; i++) {
+                     p = x * eval.V[i][k] + y * eval.V[i][k+1];
+                     if (notlast) {
+                        p = p + z * eval.V[i][k+2];
+                        eval.V[i][k+2] = eval.V[i][k+2] - p * r;
+                     }
+                     eval.V[i][k] = eval.V[i][k] - p;
+                     eval.V[i][k+1] = eval.V[i][k+1] - p * q;
+                  }
+               }  // (s != 0)
+            }  // k loop
+         }  // check convergence
+      }  // while (n >= low)
+               //if(totIter>15) std::cout<<"!!!!iter "<<totIter<<"\n";
+      
+      // Backsubstitute to find vectors of upper triangular form
+
+      if (norm == 0.0) {
+         return;
+      }
+   
+      for (n = nn-1; n >= 0; n--) {
+         p = eval.d[n];
+         q = eval.e[n];
+   
+         // float vector
+   
+         if (q == 0) {
+            int l = n;
+            eval.H[n][n] = 1.0;
+            for (int i = n-1; i >= 0; i--) {
+               w = eval.H[i][i] - p;
+               r = 0.0;
+               for (int j = l; j <= n; j++) {
+                  r = r + eval.H[i][j] * eval.H[j][n];
+               }
+               if (eval.e[i] < 0.0) {
+                  z = w;
+                  s = r;
+               } else {
+                  l = i;
+                  if (eval.e[i] == 0.0) {
+                     if (w != 0.0) {
+                        eval.H[i][n] = -r / w;
+                     } else {
+                        eval.H[i][n] = -r / (eps * norm);
+                     }
+   
+                  // Solve real equations
+   
+                  } else {
+                     x = eval.H[i][i+1];
+                     y = eval.H[i+1][i];
+                     q = (eval.d[i] - p) * (eval.d[i] - p) + eval.e[i] * eval.e[i];
+                     t = (x * s - z * r) / q;
+                     eval.H[i][n] = t;
+                     if (fabs(x) > fabs(z)) {
+                        eval.H[i+1][n] = (-r - w * t) / x;
+                     } else {
+                        eval.H[i+1][n] = (-s - y * t) / z;
+                     }
+                  }
+   
+                  // Overflow control
+   
+                  t = fabs(eval.H[i][n]);
+                  if ((eps * t) * t > 1) {
+                     for (int j = i; j <= n; j++) {
+                        eval.H[j][n] = eval.H[j][n] / t;
+                     }
+                  }
+               }
+            }
+   
+         // Complex vector
+   
+         } else if (q < 0) {
+            int l = n-1;
+
+            // Last vector component imaginary so matrix is triangular
+   
+            if (fabs(eval.H[n][n-1]) > fabs(eval.H[n-1][n])) {
+               eval.H[n-1][n-1] = q / eval.H[n][n-1];
+               eval.H[n-1][n] = -(eval.H[n][n] - p) / eval.H[n][n-1];
+            } else {
+               Eigencdiv(eval, 0.0,-eval.H[n-1][n],eval.H[n-1][n-1]-p,q);
+               eval.H[n-1][n-1] = eval.cdivr;
+               eval.H[n-1][n] = eval.cdivi;
+            }
+            eval.H[n][n-1] = 0.0;
+            eval.H[n][n] = 1.0;
+            for (int i = n-2; i >= 0; i--) {
+               float ra,sa,vr,vi;
+               ra = 0.0;
+               sa = 0.0;
+               for (int j = l; j <= n; j++) {
+                  ra = ra + eval.H[i][j] * eval.H[j][n-1];
+                  sa = sa + eval.H[i][j] * eval.H[j][n];
+               }
+               w = eval.H[i][i] - p;
+   
+               if (eval.e[i] < 0.0) {
+                  z = w;
+                  r = ra;
+                  s = sa;
+               } else {
+                  l = i;
+                  if (eval.e[i] == 0) {
+                     Eigencdiv(eval,-ra,-sa,w,q);
+                     eval.H[i][n-1] = eval.cdivr;
+                     eval.H[i][n] = eval.cdivi;
+                  } else {
+   
+                     // Solve complex equations
+   
+                     x = eval.H[i][i+1];
+                     y = eval.H[i+1][i];
+                     vr = (eval.d[i] - p) * (eval.d[i] - p) + eval.e[i] * eval.e[i] - q * q;
+                     vi = (eval.d[i] - p) * 2.0 * q;
+                     if ((vr == 0.0) && (vi == 0.0)) {
+                        vr = eps * norm * (fabs(w) + fabs(q) +
+                        fabs(x) + fabs(y) + fabs(z));
+                     }
+                     Eigencdiv(eval, x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
+                     eval.H[i][n-1] = eval.cdivr;
+                     eval.H[i][n] = eval.cdivi;
+                     if (fabs(x) > (fabs(z) + fabs(q))) {
+                        eval.H[i+1][n-1] = (-ra - w * eval.H[i][n-1] + q * eval.H[i][n]) / x;
+                        eval.H[i+1][n] = (-sa - w * eval.H[i][n] - q * eval.H[i][n-1]) / x;
+                     } else {
+                        Eigencdiv(eval, -r-y*eval.H[i][n-1],-s-y*eval.H[i][n],z,q);
+                        eval.H[i+1][n-1] = eval.cdivr;
+                        eval.H[i+1][n] = eval.cdivi;
+                     }
+                  }
+   
+                  // Overflow control
+
+                  t = max(fabs(eval.H[i][n-1]),fabs(eval.H[i][n]));
+                  if ((eps * t) * t > 1) {
+                     for (int j = i; j <= n; j++) {
+                        eval.H[j][n-1] = eval.H[j][n-1] / t;
+                        eval.H[j][n] = eval.H[j][n] / t;
+                     }
+                  }
+               }
+            }
+         }
+      }
+   
+      // Vectors of isolated roots
+   
+      for (int i = 0; i < nn; i++) {
+         if (i < low || i > high) {
+            for (int j = i; j < nn; j++) {
+               eval.V[i][j] = eval.H[i][j];
+            }
+         }
+      }
+   
+      // Back transformation to get eigenvectors of original matrix
+   
+      for (int j = nn-1; j >= low; j--) {
+         for (int i = low; i <= high; i++) {
+            z = 0.0;
+            for (int k = low; k <= min(j,high); k++) {
+               z = z + eval.V[i][k] * eval.H[k][j];
+            }
+            eval.V[i][j] = z;
+         }
+      }
+}
+
+
+
+int computeEigenvalues3x3(
+               float dout[3], 
+               float a[3][3])
+{
+  /*TNT::Array2D<float> A = TNT::Array2D<float>(3,3, &a[0][0]);
+  TNT::Array1D<float> eig = TNT::Array1D<float>(3);
+  TNT::Array1D<float> eigImag = TNT::Array1D<float>(3);
+  JAMA::Eigenvalue<float> jeig = JAMA::Eigenvalue<float>(A);*/
+
+       sEigenvalue jeig;
+
+       // Compute the values
+       {
+               jeig.n = 3;
+               int n=3;
+      //V = Array2D<float>(n,n);
+      //d = Array1D<float>(n);
+      //e = Array1D<float>(n);
+               for (int y=0; y<3; y++)
+                {
+                        jeig.d[y]=0.0f;
+                        jeig.e[y]=0.0f;
+                        for (int t=0; t<3; t++) jeig.V[y][t]=0.0f;
+                }
+
+      jeig.issymmetric = 1;
+      for (int j = 0; (j < 3) && jeig.issymmetric; j++) {
+         for (int i = 0; (i < 3) && jeig.issymmetric; i++) {
+            jeig.issymmetric = (a[i][j] == a[j][i]);
+         }
+      }
+
+      if (jeig.issymmetric) {
+         for (int i = 0; i < 3; i++) {
+            for (int j = 0; j < 3; j++) {
+               jeig.V[i][j] = a[i][j];
+            }
+         }
+   
+         // Tridiagonalize.
+         Eigentred2(jeig);
+   
+         // Diagonalize.
+         Eigentql2(jeig);
+
+      } else {
+         //H = TNT::Array2D<float>(n,n);
+            for (int y=0; y<3; y++)
+                {
+                        jeig.ort[y]=0.0f;
+                        for (int t=0; t<3; t++) jeig.H[y][t]=0.0f;
+                }
+         //ort = TNT::Array1D<float>(n);
+         
+         for (int j = 0; j < n; j++) {
+            for (int i = 0; i < n; i++) {
+               jeig.H[i][j] = a[i][j];
+            }
+         }
+   
+         // Reduce to Hessenberg form.
+         Eigenorthes(jeig);
+   
+         // Reduce Hessenberg to real Schur form.
+         Eigenhqr2(jeig);
+      }
+   }
+
+  //jeig.getfloatEigenvalues(eig);
+
+  // complex ones
+  //jeig.getImagEigenvalues(eigImag);
+  dout[0]  = sqrt(jeig.d[0]*jeig.d[0] + jeig.e[0]*jeig.e[0]);
+  dout[1]  = sqrt(jeig.d[1]*jeig.d[1] + jeig.e[1]*jeig.e[1]);
+  dout[2]  = sqrt(jeig.d[2]*jeig.d[2] + jeig.e[2]*jeig.e[2]);
+  return 0;
+}
index 6ff61c5ca8e45f344e143a5e3b1f43bfc5c017af..9c169711c09d726fbf171bc229b03060d07cb591 100644 (file)
 // along with Wavelet Turbulence.  If not, see <http://www.gnu.org/licenses/>.
 // 
 // Copyright 2008 Theodore Kim and Nils Thuerey
-// 
+//
+//////////////////////////////////////////////////////////////////////
+// Modified to not require TNT matrix library anymore. It was very slow
+// when being run in parallel. Required TNT JAMA::Eigenvalue libraries were
+// converted into independent functions.
+//             - MiikaH
+//
 //////////////////////////////////////////////////////////////////////
 // Helper function, compute eigenvalues of 3x3 matrix
 //////////////////////////////////////////////////////////////////////
 
-#include "tnt/jama_eig.h"
+#ifndef EIGENVAL_HELPER_H
+#define EIGENVAL_HELPER_H
+
+//#include "tnt/jama_eig.h"
+
+#include <algorithm>
+#include <cmath>
+
+using namespace std;
 
 //////////////////////////////////////////////////////////////////////
 // eigenvalues of 3x3 non-symmetric matrix
 //////////////////////////////////////////////////////////////////////
-int inline computeEigenvalues3x3(
-               float dout[3], 
-               float a[3][3])
+
+
+struct sEigenvalue
 {
-  TNT::Array2D<float> A = TNT::Array2D<float>(3,3, &a[0][0]);
-  TNT::Array1D<float> eig = TNT::Array1D<float>(3);
-  TNT::Array1D<float> eigImag = TNT::Array1D<float>(3);
-  JAMA::Eigenvalue<float> jeig = JAMA::Eigenvalue<float>(A);
-  jeig.getRealEigenvalues(eig);
-
-  // complex ones
-  jeig.getImagEigenvalues(eigImag);
-  dout[0]  = sqrt(eig[0]*eig[0] + eigImag[0]*eigImag[0]);
-  dout[1]  = sqrt(eig[1]*eig[1] + eigImag[1]*eigImag[1]);
-  dout[2]  = sqrt(eig[2]*eig[2] + eigImag[2]*eigImag[2]);
-  return 0;
-}
-
-#undef rfabs 
-#undef ROT
+       int n;
+       int issymmetric;
+       float d[3];         /* real part */
+       float e[3];         /* img part */
+       float V[3][3];          /* Eigenvectors */
+
+       float H[3][3];
+   
+
+    float ort[3];
+
+       float cdivr;
+       float cdivi;
+};
+
+void Eigentred2(sEigenvalue& eval);
+
+void Eigencdiv(sEigenvalue& eval, float xr, float xi, float yr, float yi);
+
+void Eigentql2 (sEigenvalue& eval);
+
+void Eigenorthes (sEigenvalue& eval);
+
+void Eigenhqr2 (sEigenvalue& eval);
+
+int computeEigenvalues3x3(float dout[3], float a[3][3]);
+
+
+#endif
index 729d73bb4f35d412f1078425014f8aa3b158e22d..25673630fc469d497243c26d5e1297f49e70d227 100644 (file)
 // FLUID_3D.cpp: implementation of the FLUID_3D class.
 //
 //////////////////////////////////////////////////////////////////////
+// Heavy parallel optimization done. Many of the old functions now
+// take begin and end parameters and process only specified part of the data.
+// Some functions were divided into multiple ones.
+//             - MiikaH
+//////////////////////////////////////////////////////////////////////
 
 #include "FLUID_3D.h"
 #include "IMAGE.h"
 #include "SPHERE.h"
 #include <zlib.h>
 
+#if PARALLEL==1
+#include <omp.h>
+#endif // PARALLEL 
+
 // boundary conditions of the fluid domain
 #define DOMAIN_BC_FRONT  0 // z
 #define DOMAIN_BC_TOP    1 // y
@@ -90,6 +99,13 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dt) :
        _heatOld      = new float[_totalCells];
        _obstacles    = new unsigned char[_totalCells]; // set 0 at end of step
 
+       // For threaded version:
+       _xVelocityTemp = new float[_totalCells];
+       _yVelocityTemp = new float[_totalCells];
+       _zVelocityTemp = new float[_totalCells];
+       _densityTemp   = new float[_totalCells];
+       _heatTemp      = new float[_totalCells];
+
        // DG TODO: check if alloc went fine
 
        for (int x = 0; x < _totalCells; x++)
@@ -167,6 +183,12 @@ FLUID_3D::~FLUID_3D()
        if (_obstacles) delete[] _obstacles;
     // if (_wTurbulence) delete _wTurbulence;
 
+       if (_xVelocityTemp) delete[] _xVelocityTemp;
+       if (_yVelocityTemp) delete[] _yVelocityTemp;
+       if (_zVelocityTemp) delete[] _zVelocityTemp;
+       if (_densityTemp) delete[] _densityTemp;
+       if (_heatTemp) delete[] _heatTemp;
+
     // printf("deleted fluid\n");
 }
 
@@ -182,108 +204,306 @@ void FLUID_3D::initBlenderRNA(float *alpha, float *beta)
 //////////////////////////////////////////////////////////////////////
 void FLUID_3D::step()
 {
-       // addSmokeTestCase(_density, _res);
-       // addSmokeTestCase(_heat, _res);
-
-       wipeBoundaries();
-
-       // run the solvers
-       addVorticity();
-       addBuoyancy(_heat, _density);
-       addForce();
-       project();
-       diffuseHeat();
-
-       // advect everything
-       advectMacCormack();
-
-       // if(_wTurbulence) {
-       //      _wTurbulence->stepTurbulenceFull(_dt/_dx,
-       //                      _xVelocity, _yVelocity, _zVelocity, _obstacles);
-               // _wTurbulence->stepTurbulenceReadable(_dt/_dx,
-               //  _xVelocity, _yVelocity, _zVelocity, _obstacles);
-       // }
-/*
- // no file output
-  float *src = _density;
-       string prefix = string("./original.preview/density_fullxy_");
-       writeImageSliceXY(src,_res, _res[2]/2, prefix, _totalSteps);
-*/
-       // artificial damping -- this is necessary because we use a
-  // collated grid, and at very coarse grid resolutions, banding
-  // artifacts can occur
-       artificialDamping(_xVelocity);
-       artificialDamping(_yVelocity);
-       artificialDamping(_zVelocity);
-/*
-// no file output
-  string pbrtPrefix = string("./pbrt/density_small_");
-  IMAGE::dumpPBRT(_totalSteps, pbrtPrefix, _density, _res[0],_res[1],_res[2]);
-  */
-       _totalTime += _dt;
-       _totalSteps++;  
 
-       // todo xxx dg: only clear obstacles, not boundaries
-       // memset(_obstacles, 0, sizeof(unsigned char)*_xRes*_yRes*_zRes);
+       int threadval = 1;
+#if PARALLEL==1
+       threadval = omp_get_max_threads();
+#endif
+
+       int stepParts = 1;
+       float partSize = _zRes;
+
+#if PARALLEL==1
+       stepParts = threadval*2;        // Dividing parallelized sections into numOfThreads * 2 sections
+       partSize = (float)_zRes/stepParts;      // Size of one part;
+
+  if (partSize < 4) {stepParts = threadval;                                    // If the slice gets too low (might actually slow things down, change it to larger
+                                       partSize = (float)_zRes/stepParts;}
+  if (partSize < 4) {stepParts = (int)(ceil((float)_zRes/4.0f));       // If it's still too low (only possible on future systems with +24 cores), change it to 4
+                                       partSize = (float)_zRes/stepParts;}
+#else
+       int zBegin=0;
+       int zEnd=_zRes;
+#endif
+
+
+#if PARALLEL==1
+       #pragma omp parallel
+       {
+       #pragma omp for schedule(static,1)
+       for (int i=0; i<stepParts; i++)
+       {
+               int zBegin = (int)((float)i*partSize + 0.5f);
+               int zEnd = (int)((float)(i+1)*partSize + 0.5f);
+#endif
+
+               wipeBoundariesSL(zBegin, zEnd);
+               addVorticity(zBegin, zEnd);
+               addBuoyancy(_heat, _density, zBegin, zEnd);
+               addForce(zBegin, zEnd);
+
+#if PARALLEL==1
+       }       // end of parallel
+       #pragma omp barrier
+
+       #pragma omp single
+       {
+#endif
+       SWAP_POINTERS(_xVelocity, _xVelocityTemp);
+       SWAP_POINTERS(_yVelocity, _yVelocityTemp);
+       SWAP_POINTERS(_zVelocity, _zVelocityTemp);
+#if PARALLEL==1
+       }       // end of single
+
+       #pragma omp barrier
+
+       #pragma omp for
+       for (int i=0; i<2; i++)
+       {
+               if (i==0)
+               {
+#endif
+               project();
+#if PARALLEL==1
+               }
+               else
+               {
+#endif
+               diffuseHeat();
+#if PARALLEL==1
+               }
+       }
+
+       #pragma omp barrier
+
+       #pragma omp single
+       {
+#endif
+       SWAP_POINTERS(_xVelocity, _xVelocityOld);
+       SWAP_POINTERS(_yVelocity, _yVelocityOld);
+       SWAP_POINTERS(_zVelocity, _zVelocityOld);
+       SWAP_POINTERS(_density, _densityOld);
+       SWAP_POINTERS(_heat, _heatOld);
+
+       advectMacCormackBegin(0, _zRes);
+
+#if PARALLEL==1
+       }       // end of single
+
+       #pragma omp barrier
+
+
+       #pragma omp for schedule(static,1)
+       for (int i=0; i<stepParts; i++)
+       {
+
+               int zBegin = (int)((float)i*partSize + 0.5f);
+               int zEnd = (int)((float)(i+1)*partSize + 0.5f);
+#endif
+
+       advectMacCormackEnd1(zBegin, zEnd);
+
+#if PARALLEL==1
+       }       // end of parallel
+
+       #pragma omp barrier
+
+       #pragma omp for schedule(static,1)
+       for (int i=0; i<stepParts; i++)
+       {
+
+               int zBegin = (int)((float)i*partSize + 0.5f);
+               int zEnd = (int)((float)(i+1)*partSize + 0.5f);
+#endif
+
+               advectMacCormackEnd2(zBegin, zEnd);
+
+               artificialDampingSL(zBegin, zEnd);
+
+               // Using forces as temp arrays
+
+#if PARALLEL==1
+       }
+       }
+
+
+
+       for (int i=1; i<stepParts; i++)
+       {
+               int zPos=(int)((float)i*partSize + 0.5f);
+               
+               artificialDampingExactSL(zPos);
+
+       }
+#endif
+
+       SWAP_POINTERS(_xVelocity, _xForce);
+       SWAP_POINTERS(_yVelocity, _yForce);
+       SWAP_POINTERS(_zVelocity, _zForce);
+
+
+
+
+       _totalTime += _dt;
+       _totalSteps++;          
 
-       // wipe forces
-       // for external forces we can't do it at the beginning of this function but at the end
        for (int i = 0; i < _totalCells; i++)
        {
                _xForce[i] = _yForce[i] = _zForce[i] = 0.0f;
        }
+
 }
 
 //////////////////////////////////////////////////////////////////////
 // helper function to dampen co-located grid artifacts of given arrays in intervals
 // (only needed for velocity, strength (w) depends on testcase...
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::artificialDamping(float* field) {
+
+
+void FLUID_3D::artificialDampingSL(int zBegin, int zEnd) {
        const float w = 0.9;
+
+       memmove(_xForce+(_slabSize*zBegin), _xVelocityTemp+(_slabSize*zBegin), sizeof(float)*_slabSize*(zEnd-zBegin));
+       memmove(_yForce+(_slabSize*zBegin), _yVelocityTemp+(_slabSize*zBegin), sizeof(float)*_slabSize*(zEnd-zBegin));
+       memmove(_zForce+(_slabSize*zBegin), _zVelocityTemp+(_slabSize*zBegin), sizeof(float)*_slabSize*(zEnd-zBegin));
+
+
        if(_totalSteps % 4 == 1) {
-               for (int z = 1; z < _res[2]-1; z++)
+               for (int z = zBegin+1; z < zEnd-1; z++)
                        for (int y = 1; y < _res[1]-1; y++)
                                for (int x = 1+(y+z)%2; x < _res[0]-1; x+=2) {
                                        const int index = x + y*_res[0] + z * _slabSize;
-                                       field[index] = (1-w)*field[index] + 1./6. * w*(
-                                                       field[index+1] + field[index-1] +
-                                                       field[index+_res[0]] + field[index-_res[0]] +
-                                                       field[index+_slabSize] + field[index-_slabSize] );
+                                       _xForce[index] = (1-w)*_xVelocityTemp[index] + 1./6. * w*(
+                                                       _xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
+                                                       _xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
+                                                       _xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
+
+                                       _yForce[index] = (1-w)*_yVelocityTemp[index] + 1./6. * w*(
+                                                       _yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
+                                                       _yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
+                                                       _yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
+
+                                       _zForce[index] = (1-w)*_zVelocityTemp[index] + 1./6. * w*(
+                                                       _zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
+                                                       _zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
+                                                       _zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
                                }
        }
+
        if(_totalSteps % 4 == 3) {
-               for (int z = 1; z < _res[2]-1; z++)
+               for (int z = zBegin+1; z < zEnd-1; z++)
                        for (int y = 1; y < _res[1]-1; y++)
                                for (int x = 1+(y+z+1)%2; x < _res[0]-1; x+=2) {
                                        const int index = x + y*_res[0] + z * _slabSize;
-                                       field[index] = (1-w)*field[index] + 1./6. * w*(
-                                                       field[index+1] + field[index-1] +
-                                                       field[index+_res[0]] + field[index-_res[0]] +
-                                                       field[index+_slabSize] + field[index-_slabSize] );
+                                       _xForce[index] = (1-w)*_xVelocityTemp[index] + 1./6. * w*(
+                                                       _xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
+                                                       _xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
+                                                       _xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
+
+                                       _yForce[index] = (1-w)*_yVelocityTemp[index] + 1./6. * w*(
+                                                       _yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
+                                                       _yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
+                                                       _yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
+
+                                       _zForce[index] = (1-w)*_zVelocityTemp[index] + 1./6. * w*(
+                                                       _zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
+                                                       _zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
+                                                       _zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
                                }
+
+       }
+}
+
+
+
+void FLUID_3D::artificialDampingExactSL(int pos) {
+       const float w = 0.9;
+       int index, x,y,z;
+       
+
+       size_t posslab;
+
+       for (z=pos-1; z<=pos; z++)
+       {
+       posslab=z * _slabSize;
+
+       if(_totalSteps % 4 == 1) {
+                       for (y = 1; y < _res[1]-1; y++)
+                               for (x = 1+(y+z)%2; x < _res[0]-1; x+=2) {
+                                       index = x + y*_res[0] + posslab;
+                                       _xForce[index] = (1-w)*_xVelocityTemp[index] + 1./6. * w*(
+                                                       _xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
+                                                       _xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
+                                                       _xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
+
+                                       _yForce[index] = (1-w)*_yVelocityTemp[index] + 1./6. * w*(
+                                                       _yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
+                                                       _yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
+                                                       _yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
+
+                                       _zForce[index] = (1-w)*_zVelocityTemp[index] + 1./6. * w*(
+                                                       _zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
+                                                       _zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
+                                                       _zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
+                                       
+                               }
+       }
+
+       if(_totalSteps % 4 == 3) {
+                       for (y = 1; y < _res[1]-1; y++)
+                               for (x = 1+(y+z+1)%2; x < _res[0]-1; x+=2) {
+                                       index = x + y*_res[0] + posslab;
+                                       _xForce[index] = (1-w)*_xVelocityTemp[index] + 1./6. * w*(
+                                                       _xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
+                                                       _xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
+                                                       _xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
+
+                                       _yForce[index] = (1-w)*_yVelocityTemp[index] + 1./6. * w*(
+                                                       _yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
+                                                       _yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
+                                                       _yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
+
+                                       _zForce[index] = (1-w)*_zVelocityTemp[index] + 1./6. * w*(
+                                                       _zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
+                                                       _zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
+                                                       _zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
+                                       
+                               }
+
+       }
        }
 }
 
 //////////////////////////////////////////////////////////////////////
 // copy out the boundary in all directions
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::copyBorderAll(float* field)
+void FLUID_3D::copyBorderAll(float* field, int zBegin, int zEnd)
 {
-       int index;
+       int index, x, y, z;
+       int zSize = zEnd-zBegin;
+       int _blockTotalCells=_slabSize * zSize;
+
+       if ((zBegin==0))
        for (int y = 0; y < _yRes; y++)
                for (int x = 0; x < _xRes; x++)
                {
                        // front slab
                        index = x + y * _xRes;
                        field[index] = field[index + _slabSize];
+    }
+
+       if (zEnd==_zRes)
+       for (y = 0; y < _yRes; y++)
+               for (x = 0; x < _xRes; x++)
+               {
 
                        // back slab
-                       index += _totalCells - _slabSize;
+                       index = x + y * _xRes + _blockTotalCells - _slabSize;
                        field[index] = field[index - _slabSize];
     }
 
-       for (int z = 0; z < _zRes; z++)
-               for (int x = 0; x < _xRes; x++)
+       for (z = 0; z < zSize; z++)
+               for (x = 0; x < _xRes; x++)
     {
                        // bottom slab
                        index = x + z * _slabSize;
@@ -294,8 +514,8 @@ void FLUID_3D::copyBorderAll(float* field)
                        field[index] = field[index - _xRes];
     }
 
-       for (int z = 0; z < _zRes; z++)
-               for (int y = 0; y < _yRes; y++)
+       for (z = 0; z < zSize; z++)
+               for (y = 0; y < _yRes; y++)
     {
                        // left slab
                        index = y * _xRes + z * _slabSize;
@@ -310,27 +530,123 @@ void FLUID_3D::copyBorderAll(float* field)
 //////////////////////////////////////////////////////////////////////
 // wipe boundaries of velocity and density
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::wipeBoundaries()
+void FLUID_3D::wipeBoundaries(int zBegin, int zEnd)
 {
-       setZeroBorder(_xVelocity, _res);
-       setZeroBorder(_yVelocity, _res);
-       setZeroBorder(_zVelocity, _res);
-       setZeroBorder(_density, _res);
+       setZeroBorder(_xVelocity, _res, zBegin, zEnd);
+       setZeroBorder(_yVelocity, _res, zBegin, zEnd);
+       setZeroBorder(_zVelocity, _res, zBegin, zEnd);
+       setZeroBorder(_density, _res, zBegin, zEnd);
 }
 
+void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd)
+{
+       
+       /////////////////////////////////////
+       // setZeroBorder to all:
+       /////////////////////////////////////
+
+       /////////////////////////////////////
+       // setZeroX
+       /////////////////////////////////////
+
+       const int slabSize = _xRes * _yRes;
+       int index, x,y,z;
+
+       for (z = zBegin; z < zEnd; z++)
+               for (y = 0; y < _yRes; y++)
+               {
+                       // left slab
+                       index = y * _xRes + z * slabSize;
+                       _xVelocity[index] = 0.0f;
+                       _yVelocity[index] = 0.0f;
+                       _zVelocity[index] = 0.0f;
+                       _density[index] = 0.0f;
+
+                       // right slab
+                       index += _xRes - 1;
+                       _xVelocity[index] = 0.0f;
+                       _yVelocity[index] = 0.0f;
+                       _zVelocity[index] = 0.0f;
+                       _density[index] = 0.0f;
+               }
+
+       /////////////////////////////////////
+       // setZeroY
+       /////////////////////////////////////
+
+       for (z = zBegin; z < zEnd; z++)
+               for (x = 0; x < _xRes; x++)
+               {
+                       // bottom slab
+                       index = x + z * slabSize;
+                       _xVelocity[index] = 0.0f;
+                       _yVelocity[index] = 0.0f;
+                       _zVelocity[index] = 0.0f;
+                       _density[index] = 0.0f;
+
+                       // top slab
+                       index += slabSize - _xRes;
+                       _xVelocity[index] = 0.0f;
+                       _yVelocity[index] = 0.0f;
+                       _zVelocity[index] = 0.0f;
+                       _density[index] = 0.0f;
+
+               }
+
+       /////////////////////////////////////
+       // setZeroZ
+       /////////////////////////////////////
+
+
+       const int totalCells = _xRes * _yRes * _zRes;
+
+       index = 0;
+       if ((zBegin == 0))
+       for (y = 0; y < _yRes; y++)
+               for (x = 0; x < _xRes; x++, index++)
+               {
+                       // front slab
+                       _xVelocity[index] = 0.0f;
+                       _yVelocity[index] = 0.0f;
+                       _zVelocity[index] = 0.0f;
+                       _density[index] = 0.0f;
+    }
+
+       if (zEnd == _zRes)
+       {
+               index=0;
+               int indexx=0;
+               const int cellsslab = totalCells - slabSize;
+
+               for (y = 0; y < _yRes; y++)
+                       for (x = 0; x < _xRes; x++, index++)
+                       {
+
+                               // back slab
+                               indexx = index + cellsslab;
+                               _xVelocity[indexx] = 0.0f;
+                               _yVelocity[indexx] = 0.0f;
+                               _zVelocity[indexx] = 0.0f;
+                               _density[indexx] = 0.0f;
+                       }
+       }
+
+}
 //////////////////////////////////////////////////////////////////////
 // add forces to velocity field
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::addForce()
+void FLUID_3D::addForce(int zBegin, int zEnd)
 {
-       for (int i = 0; i < _totalCells; i++)
+       int begin=zBegin * _slabSize;
+       int end=begin + (zEnd - zBegin) * _slabSize;
+
+       for (int i = begin; i < end; i++)
        {
-               _xVelocity[i] += _dt * _xForce[i];
-               _yVelocity[i] += _dt * _yForce[i];
-               _zVelocity[i] += _dt * _zForce[i];
+               _xVelocityTemp[i] = _xVelocity[i] + _dt * _xForce[i];
+               _yVelocityTemp[i] = _yVelocity[i] + _dt * _yForce[i];
+               _zVelocityTemp[i] = _zVelocity[i] + _dt * _zForce[i];
        }
 }
-
 //////////////////////////////////////////////////////////////////////
 // project into divergence free field
 //////////////////////////////////////////////////////////////////////
@@ -344,18 +660,18 @@ void FLUID_3D::project()
 
        memset(_pressure, 0, sizeof(float)*_totalCells);
        memset(_divergence, 0, sizeof(float)*_totalCells);
-
-       setObstacleBoundaries(_pressure);
-
+       
+       setObstacleBoundaries(_pressure, 0, _zRes);
+       
        // copy out the boundaries
-       if(DOMAIN_BC_LEFT == 0)  setNeumannX(_xVelocity, _res);
-       else setZeroX(_xVelocity, _res); 
+       if(DOMAIN_BC_LEFT == 0)  setNeumannX(_xVelocity, _res, 0, _zRes);
+       else setZeroX(_xVelocity, _res, 0, _zRes); 
 
-       if(DOMAIN_BC_TOP == 0)   setNeumannY(_yVelocity, _res);
-       else setZeroY(_yVelocity, _res); 
+       if(DOMAIN_BC_TOP == 0)   setNeumannY(_yVelocity, _res, 0, _zRes);
+       else setZeroY(_yVelocity, _res, 0, _zRes); 
 
-       if(DOMAIN_BC_FRONT == 0) setNeumannZ(_zVelocity, _res);
-       else setZeroZ(_zVelocity, _res);
+       if(DOMAIN_BC_FRONT == 0) setNeumannZ(_zVelocity, _res, 0, _zRes);
+       else setZeroZ(_zVelocity, _res, 0, _zRes);
 
        // calculate divergence
        index = _slabSize + _xRes + 1;
@@ -385,12 +701,12 @@ void FLUID_3D::project()
                                // DG: commenting this helps CG to get a better start, 10-20% speed improvement
                                // _pressure[index] = 0.0f;
                        }
-       copyBorderAll(_pressure);
+       copyBorderAll(_pressure, 0, _zRes);
 
        // solve Poisson equation
        solvePressurePre(_pressure, _divergence, _obstacles);
 
-       setObstaclePressure(_pressure);
+       setObstaclePressure(_pressure, 0, _zRes);
 
        // project out solution
        float invDx = 1.0f / _dx;
@@ -411,6 +727,9 @@ void FLUID_3D::project()
        if (_divergence) delete[] _divergence;
 }
 
+
+
+
 //////////////////////////////////////////////////////////////////////
 // diffuse heat
 //////////////////////////////////////////////////////////////////////
@@ -418,13 +737,14 @@ void FLUID_3D::diffuseHeat()
 {
        SWAP_POINTERS(_heat, _heatOld);
 
-       copyBorderAll(_heatOld);
+       copyBorderAll(_heatOld, 0, _zRes);
        solveHeat(_heat, _heatOld, _obstacles);
 
        // zero out inside obstacles
        for (int x = 0; x < _totalCells; x++)
                if (_obstacles[x])
                        _heat[x] = 0.0f;
+
 }
 
 //////////////////////////////////////////////////////////////////////
@@ -444,12 +764,28 @@ void FLUID_3D::addObstacle(OBSTACLE* obstacle)
 //////////////////////////////////////////////////////////////////////
 // calculate the obstacle directional types
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::setObstaclePressure(float *_pressure)
+void FLUID_3D::setObstaclePressure(float *_pressure, int zBegin, int zEnd)
 {
+
+       // compleately TODO
+
+       const size_t index_ = _slabSize + _xRes + 1;
+
+       //int vIndex=_slabSize + _xRes + 1;
+
+       int bb=0;
+       int bt=0;
+
+       if (zBegin == 0) {bb = 1;}
+       if (zEnd == _zRes) {bt = 1;}
+
        // tag remaining obstacle blocks
-       for (int z = 1, index = _slabSize + _xRes + 1;
-                       z < _zRes - 1; z++, index += 2 * _xRes)
+       for (int z = zBegin + bb; z < zEnd - bt; z++)
+       {
+               size_t index = index_ +(z-1)*_slabSize;
+
                for (int y = 1; y < _yRes - 1; y++, index += 2)
+               {
                        for (int x = 1; x < _xRes - 1; x++, index++)
                {
                        // could do cascade of ifs, but they are a pain
@@ -507,15 +843,33 @@ void FLUID_3D::setObstaclePressure(float *_pressure)
                                // this means it's not a full no-slip boundary condition
                                // but a "half-slip" - still looks ok right now
                        }
-               }
+                       //vIndex++;
+               }       // x loop
+               //vIndex += 2;
+               }       // y loop
+               //vIndex += 2 * _xRes;
+       }       // z loop
 }
 
-void FLUID_3D::setObstacleBoundaries(float *_pressure)
+void FLUID_3D::setObstacleBoundaries(float *_pressure, int zBegin, int zEnd)
 {
        // cull degenerate obstacles , move to addObstacle?
-       for (int z = 1, index = _slabSize + _xRes + 1;
-                       z < _zRes - 1; z++, index += 2 * _xRes)
+
+       // r = b - Ax
+       const size_t index_ = _slabSize + _xRes + 1;
+
+       int bb=0;
+       int bt=0;
+
+       if (zBegin == 0) {bb = 1;}
+       if (zEnd == _zRes) {bt = 1;}
+
+       for (int z = zBegin + bb; z < zEnd - bt; z++)
+       {
+               size_t index = index_ +(z-1)*_slabSize;
+               
                for (int y = 1; y < _yRes - 1; y++, index += 2)
+               {
                        for (int x = 1; x < _xRes - 1; x++, index++)
                        {
                                if (_obstacles[index] != EMPTY)
@@ -545,17 +899,22 @@ void FLUID_3D::setObstacleBoundaries(float *_pressure)
                                        _zVelocity[index] = 0.0f;
                                        _pressure[index] = 0.0f;
                                }
-                       }
+                               //vIndex++;
+                       }       // x-loop
+                       //vIndex += 2;
+               }       // y-loop
+               //vIndex += 2* _xRes;
+       }       // z-loop
 }
 
 //////////////////////////////////////////////////////////////////////
 // add buoyancy forces
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::addBuoyancy(float *heat, float *density)
+void FLUID_3D::addBuoyancy(float *heat, float *density, int zBegin, int zEnd)
 {
-       int index = 0;
+       int index = zBegin*_slabSize;
 
-       for (int z = 0; z < _zRes; z++)
+       for (int z = zBegin; z < zEnd; z++)
                for (int y = 0; y < _yRes; y++)
                        for (int x = 0; x < _xRes; x++, index++)
                        {
@@ -566,30 +925,55 @@ void FLUID_3D::addBuoyancy(float *heat, float *density)
 //////////////////////////////////////////////////////////////////////
 // add vorticity to the force field
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::addVorticity()
+void FLUID_3D::addVorticity(int zBegin, int zEnd)
 {
-       int x,y,z,index;
+       //int x,y,z,index;
        if(_vorticityEps<=0.) return;
 
+       int _blockSize=zEnd-zBegin;
+       int _blockTotalCells = _slabSize * (_blockSize+2);
+
        float *_xVorticity, *_yVorticity, *_zVorticity, *_vorticity;
 
-       _xVorticity = new float[_totalCells];
-       _yVorticity = new float[_totalCells];
-       _zVorticity = new float[_totalCells];
-       _vorticity = new float[_totalCells];
+       int _vIndex = _slabSize + _xRes + 1;
+       int bb=0;
+       int bt=0;
+       int bb1=-1;
+       int bt1=-1;
+
+       if (zBegin == 0) {bb1 = 1; bb = 1; _blockTotalCells-=_blockSize;}
+       if (zEnd == _zRes) {bt1 = 1;bt = 1; _blockTotalCells-=_blockSize;}
+
+       _xVorticity = new float[_blockTotalCells];
+       _yVorticity = new float[_blockTotalCells];
+       _zVorticity = new float[_blockTotalCells];
+       _vorticity = new float[_blockTotalCells];
 
-       memset(_xVorticity, 0, sizeof(float)*_totalCells);
-       memset(_yVorticity, 0, sizeof(float)*_totalCells);
-       memset(_zVorticity, 0, sizeof(float)*_totalCells);
-       memset(_vorticity, 0, sizeof(float)*_totalCells);
+       memset(_xVorticity, 0, sizeof(float)*_blockTotalCells);
+       memset(_yVorticity, 0, sizeof(float)*_blockTotalCells);
+       memset(_zVorticity, 0, sizeof(float)*_blockTotalCells);
+       memset(_vorticity, 0, sizeof(float)*_blockTotalCells);
+
+       //const size_t indexsetupV=_slabSize;
+       const size_t index_ = _slabSize + _xRes + 1;
 
        // calculate vorticity
        float gridSize = 0.5f / _dx;
-       index = _slabSize + _xRes + 1;
-       for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-               for (y = 1; y < _yRes - 1; y++, index += 2)
-                       for (x = 1; x < _xRes - 1; x++, index++)
+       //index = _slabSize + _xRes + 1;
+
+
+       size_t vIndex=_xRes + 1;
+
+       for (int z = zBegin + bb1; z < (zEnd - bt1); z++)
+       {
+               size_t index = index_ +(z-1)*_slabSize;
+               vIndex = index-(zBegin-1+bb)*_slabSize;
+
+               for (int y = 1; y < _yRes - 1; y++, index += 2)
+               {
+                       for (int x = 1; x < _xRes - 1; x++, index++)
                        {
+
                                int up    = _obstacles[index + _xRes] ? index : index + _xRes;
                                int down  = _obstacles[index - _xRes] ? index : index - _xRes;
                                float dy  = (up == index || down == index) ? 1.0f / _dx : gridSize;
@@ -600,34 +984,51 @@ void FLUID_3D::addVorticity()
                                int left  = _obstacles[index - 1] ? index : index - 1;
                                float dx  = (right == index || right == index) ? 1.0f / _dx : gridSize;
 
-                               _xVorticity[index] = (_zVelocity[up] - _zVelocity[down]) * dy + (-_yVelocity[out] + _yVelocity[in]) * dz;
-                               _yVorticity[index] = (_xVelocity[out] - _xVelocity[in]) * dz + (-_zVelocity[right] + _zVelocity[left]) * dx;
-                               _zVorticity[index] = (_yVelocity[right] - _yVelocity[left]) * dx + (-_xVelocity[up] + _xVelocity[down])* dy;
+                               _xVorticity[vIndex] = (_zVelocity[up] - _zVelocity[down]) * dy + (-_yVelocity[out] + _yVelocity[in]) * dz;
+                               _yVorticity[vIndex] = (_xVelocity[out] - _xVelocity[in]) * dz + (-_zVelocity[right] + _zVelocity[left]) * dx;
+                               _zVorticity[vIndex] = (_yVelocity[right] - _yVelocity[left]) * dx + (-_xVelocity[up] + _xVelocity[down])* dy;
 
-                               _vorticity[index] = sqrtf(_xVorticity[index] * _xVorticity[index] +
-                                               _yVorticity[index] * _yVorticity[index] +
-                                               _zVorticity[index] * _zVorticity[index]);
+                               _vorticity[vIndex] = sqrtf(_xVorticity[vIndex] * _xVorticity[vIndex] +
+                                               _yVorticity[vIndex] * _yVorticity[vIndex] +
+                                               _zVorticity[vIndex] * _zVorticity[vIndex]);
+
+                               vIndex++;
                        }
+                       vIndex+=2;
+               }
+               //vIndex+=2*_xRes;
+       }
 
        // calculate normalized vorticity vectors
        float eps = _vorticityEps;
-       index = _slabSize + _xRes + 1;
-       for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-               for (y = 1; y < _yRes - 1; y++, index += 2)
-                       for (x = 1; x < _xRes - 1; x++, index++)
+       
+       //index = _slabSize + _xRes + 1;
+       vIndex=_slabSize + _xRes + 1;
+
+       for (int z = zBegin + bb; z < (zEnd - bt); z++)
+       {
+               size_t index = index_ +(z-1)*_slabSize;
+               vIndex = index-(zBegin-1+bb)*_slabSize;
+
+               for (int y = 1; y < _yRes - 1; y++, index += 2)
+               {
+                       for (int x = 1; x < _xRes - 1; x++, index++)
+                       {
+                               //
+
                                if (!_obstacles[index])
                                {
                                        float N[3];
 
-                                       int up    = _obstacles[index + _xRes] ? index : index + _xRes;
-                                       int down  = _obstacles[index - _xRes] ? index : index - _xRes;
-                                       float dy  = (up == index || down == index) ? 1.0f / _dx : gridSize;
-                                       int out   = _obstacles[index + _slabSize] ? index : index + _slabSize;
-                                       int in    = _obstacles[index - _slabSize] ? index : index - _slabSize;
-                                       float dz  = (out == index || in == index) ? 1.0f / _dx : gridSize;
-                                       int right = _obstacles[index + 1] ? index : index + 1;
-                                       int left  = _obstacles[index - 1] ? index : index - 1;
-                                       float dx  = (right == index || right == index) ? 1.0f / _dx : gridSize;
+                                       int up    = _obstacles[index + _xRes] ? vIndex : vIndex + _xRes;
+                                       int down  = _obstacles[index - _xRes] ? vIndex : vIndex - _xRes;
+                                       float dy  = (up == vIndex || down == vIndex) ? 1.0f / _dx : gridSize;
+                                       int out   = _obstacles[index + _slabSize] ? vIndex : vIndex + _slabSize;
+                                       int in    = _obstacles[index - _slabSize] ? vIndex : vIndex - _slabSize;
+                                       float dz  = (out == vIndex || in == vIndex) ? 1.0f / _dx : gridSize;
+                                       int right = _obstacles[index + 1] ? vIndex : vIndex + 1;
+                                       int left  = _obstacles[index - 1] ? vIndex : vIndex - 1;
+                                       float dx  = (right == vIndex || right == vIndex) ? 1.0f / _dx : gridSize;
                                        N[0] = (_vorticity[right] - _vorticity[left]) * dx;
                                        N[1] = (_vorticity[up] - _vorticity[down]) * dy;
                                        N[2] = (_vorticity[out] - _vorticity[in]) * dz;
@@ -640,11 +1041,17 @@ void FLUID_3D::addVorticity()
                                                N[1] *= magnitude;
                                                N[2] *= magnitude;
 
-                                               _xForce[index] += (N[1] * _zVorticity[index] - N[2] * _yVorticity[index]) * _dx * eps;
-                                               _yForce[index] -= (N[0] * _zVorticity[index] - N[2] * _xVorticity[index]) * _dx * eps;
-                                               _zForce[index] += (N[0] * _yVorticity[index] - N[1] * _xVorticity[index]) * _dx * eps;
+                                               _xForce[index] += (N[1] * _zVorticity[vIndex] - N[2] * _yVorticity[vIndex]) * _dx * eps;
+                                               _yForce[index] -= (N[0] * _zVorticity[vIndex] - N[2] * _xVorticity[vIndex]) * _dx * eps;
+                                               _zForce[index] += (N[0] * _yVorticity[vIndex] - N[1] * _xVorticity[vIndex]) * _dx * eps;
                                        }
-                               }
+                                       }       // if
+                                       vIndex++;
+                                       }       // x loop
+                               vIndex+=2;
+                               }               // y loop
+                       //vIndex+=2*_xRes;
+               }                               // z loop
                                
        if (_xVorticity) delete[] _xVorticity;
        if (_yVorticity) delete[] _yVorticity;
@@ -652,54 +1059,80 @@ void FLUID_3D::addVorticity()
        if (_vorticity) delete[] _vorticity;
 }
 
+
+void FLUID_3D::advectMacCormackBegin(int zBegin, int zEnd)
+{
+       Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
+
+       if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocityOld, res, zBegin, zEnd);
+       else setZeroX(_xVelocityOld, res, zBegin, zEnd);
+
+       if(DOMAIN_BC_TOP == 0) copyBorderY(_yVelocityOld, res, zBegin, zEnd);
+       else setZeroY(_yVelocityOld, res, zBegin, zEnd); 
+
+       if(DOMAIN_BC_FRONT == 0) copyBorderZ(_zVelocityOld, res, zBegin, zEnd);
+       else setZeroZ(_zVelocityOld, res, zBegin, zEnd);
+}
+
 //////////////////////////////////////////////////////////////////////
 // Advect using the MacCormack method from the Selle paper
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::advectMacCormack()
+void FLUID_3D::advectMacCormackEnd1(int zBegin, int zEnd)
 {
        Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
 
-       if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocity, res);
-       else setZeroX(_xVelocity, res); 
+       const float dt0 = _dt / _dx;
 
-       if(DOMAIN_BC_TOP == 0) copyBorderY(_yVelocity, res);
-       else setZeroY(_yVelocity, res); 
+       int begin=zBegin * _slabSize;
+       int end=begin + (zEnd - zBegin) * _slabSize;
+       for (int x = begin; x < end; x++)
+               _xForce[x] = 0.0;
 
-       if(DOMAIN_BC_FRONT == 0) copyBorderZ(_zVelocity, res);
-       else setZeroZ(_zVelocity, res);
+       // advectFieldMacCormack1(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res)
 
-       SWAP_POINTERS(_xVelocity, _xVelocityOld);
-       SWAP_POINTERS(_yVelocity, _yVelocityOld);
-       SWAP_POINTERS(_zVelocity, _zVelocityOld);
-       SWAP_POINTERS(_density, _densityOld);
-       SWAP_POINTERS(_heat, _heatOld);
+       advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _densityTemp, res, zBegin, zEnd);
+       advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heatTemp, res, zBegin, zEnd);
+       advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _xVelocityOld, _xVelocity, res, zBegin, zEnd);
+       advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _yVelocityOld, _yVelocity, res, zBegin, zEnd);
+       advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocity, res, zBegin, zEnd);
 
+       // Have to wait untill all the threads are done -> so continuing in step 3
+}
+
+//////////////////////////////////////////////////////////////////////
+// Advect using the MacCormack method from the Selle paper
+//////////////////////////////////////////////////////////////////////
+void FLUID_3D::advectMacCormackEnd2(int zBegin, int zEnd)
+{
        const float dt0 = _dt / _dx;
-       // use force arrays as temp arrays
-       for (int x = 0; x < _totalCells; x++)
-               _xForce[x] = _yForce[x] = 0.0;
+       Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
 
+       // use force array as temp arrays
        float* t1 = _xForce;
-       float* t2 = _yForce;
 
-       advectFieldMacCormack(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _density, t1,t2, res, _obstacles);
-       advectFieldMacCormack(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heat, t1,t2, res, _obstacles);
-       advectFieldMacCormack(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _xVelocityOld, _xVelocity, t1,t2, res, _obstacles);
-       advectFieldMacCormack(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _yVelocityOld, _yVelocity, t1,t2, res, _obstacles);
-       advectFieldMacCormack(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocity, t1,t2, res, _obstacles);
+       // advectFieldMacCormack2(dt, xVelocity, yVelocity, zVelocity, oldField, newField, tempfield, temp, res, obstacles)
+
+       advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _density, _densityTemp, t1, res, _obstacles, zBegin, zEnd);
+       advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heat, _heatTemp, t1, res, _obstacles, zBegin, zEnd);
+       advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _xVelocityOld, _xVelocityTemp, _xVelocity, t1, res, _obstacles, zBegin, zEnd);
+       advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _yVelocityOld, _yVelocityTemp, _yVelocity, t1, res, _obstacles, zBegin, zEnd);
+       advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocityTemp, _zVelocity, t1, res, _obstacles, zBegin, zEnd);
+
+       if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocityTemp, res, zBegin, zEnd);
+       else setZeroX(_xVelocityTemp, res, zBegin, zEnd);                               
 
-       if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocity, res);
-       else setZeroX(_xVelocity, res); 
+       if(DOMAIN_BC_TOP == 0) copyBorderY(_yVelocityTemp, res, zBegin, zEnd);
+       else setZeroY(_yVelocityTemp, res, zBegin, zEnd); 
 
-       if(DOMAIN_BC_TOP == 0) copyBorderY(_yVelocity, res);
-       else setZeroY(_yVelocity, res); 
+       if(DOMAIN_BC_FRONT == 0) copyBorderZ(_zVelocityTemp, res, zBegin, zEnd);
+       else setZeroZ(_zVelocityTemp, res, zBegin, zEnd);
 
-       if(DOMAIN_BC_FRONT == 0) copyBorderZ(_zVelocity, res);
-       else setZeroZ(_zVelocity, res);
+       setZeroBorder(_density, res, zBegin, zEnd);
+       setZeroBorder(_heat, res, zBegin, zEnd);
 
-       setZeroBorder(_density, res);
-       setZeroBorder(_heat, res);
 
-  for (int x = 0; x < _totalCells; x++)
-    t1[x] = t2[x] = 0.0;
+       /*int begin=zBegin * _slabSize;
+       int end=begin + (zEnd - zBegin) * _slabSize;
+  for (int x = begin; x < end; x++)
+    _xForce[x] = _yForce[x] = 0.0f;*/
 }
index a7be7f58335d38b4ca5c1704f706df96862b3361..ab121ef506eff24a1f5f328b99aae9a31b462051 100644 (file)
 // FLUID_3D.h: interface for the FLUID_3D class.
 //
 //////////////////////////////////////////////////////////////////////
+// Heavy parallel optimization done. Many of the old functions now
+// take begin and end parameters and process only specified part of the data.
+// Some functions were divided into multiple ones.
+//             - MiikaH
+//////////////////////////////////////////////////////////////////////
 
 #ifndef FLUID_3D_H
 #define FLUID_3D_H
@@ -74,7 +79,8 @@ class FLUID_3D
                int _totalImgDumps;
                int _totalVelDumps;
 
-    void artificialDamping(float* field);
+               void artificialDampingSL(int zBegin, int zEnd);
+               void artificialDampingExactSL(int pos);
 
                // fields
                float* _density;
@@ -92,6 +98,13 @@ class FLUID_3D
                float* _zForce;
                unsigned char*  _obstacles;
 
+               // Required for proper threading:
+               float* _xVelocityTemp;
+               float* _yVelocityTemp;
+               float* _zVelocityTemp;
+               float* _heatTemp;
+               float* _densityTemp;
+
                // CG fields
                int _iterations;
 
@@ -107,13 +120,14 @@ class FLUID_3D
                // WTURBULENCE* _wTurbulence;
 
                // boundary setting functions
-               void copyBorderAll(float* field);
+               void copyBorderAll(float* field, int zBegin, int zEnd);
 
                // timestepping functions
-               void wipeBoundaries();
-               void addForce();
-               void addVorticity();
-               void addBuoyancy(float *heat, float *density);
+               void wipeBoundaries(int zBegin, int zEnd);
+               void wipeBoundariesSL(int zBegin, int zEnd);
+               void addForce(int zBegin, int zEnd);
+               void addVorticity(int zBegin, int zEnd);
+               void addBuoyancy(float *heat, float *density, int zBegin, int zEnd);
 
                // solver stuff
                void project();
@@ -122,41 +136,58 @@ class FLUID_3D
                void solvePressurePre(float* field, float* b, unsigned char* skip);
                void solveHeat(float* field, float* b, unsigned char* skip);
 
+
                // handle obstacle boundaries
-               void setObstacleBoundaries(float *_pressure);
-               void setObstaclePressure(float *_pressure);
+               void setObstacleBoundaries(float *_pressure, int zBegin, int zEnd);
+               void setObstaclePressure(float *_pressure, int zBegin, int zEnd);
 
        public:
                // advection, accessed e.g. by WTURBULENCE class
-               void advectMacCormack();
+               //void advectMacCormack();
+               void advectMacCormackBegin(int zBegin, int zEnd);
+               void advectMacCormackEnd1(int zBegin, int zEnd);
+               void advectMacCormackEnd2(int zBegin, int zEnd);
 
                // boundary setting functions
-               static void copyBorderX(float* field, Vec3Int res);
-               static void copyBorderY(float* field, Vec3Int res);
-               static void copyBorderZ(float* field, Vec3Int res);
-               static void setNeumannX(float* field, Vec3Int res);
-               static void setNeumannY(float* field, Vec3Int res);
-               static void setNeumannZ(float* field, Vec3Int res);
-               static void setZeroX(float* field, Vec3Int res);
-               static void setZeroY(float* field, Vec3Int res);
-               static void setZeroZ(float* field, Vec3Int res);
-               static void setZeroBorder(float* field, Vec3Int res) {
-                       setZeroX(field, res);
-                       setZeroY(field, res);
-                       setZeroZ(field, res);
+               static void copyBorderX(float* field, Vec3Int res, int zBegin, int zEnd);
+               static void copyBorderY(float* field, Vec3Int res, int zBegin, int zEnd);
+               static void copyBorderZ(float* field, Vec3Int res, int zBegin, int zEnd);
+               static void setNeumannX(float* field, Vec3Int res, int zBegin, int zEnd);
+               static void setNeumannY(float* field, Vec3Int res, int zBegin, int zEnd);
+               static void setNeumannZ(float* field, Vec3Int res, int zBegin, int zEnd);
+               static void setZeroX(float* field, Vec3Int res, int zBegin, int zEnd);
+               static void setZeroY(float* field, Vec3Int res, int zBegin, int zEnd);
+               static void setZeroZ(float* field, Vec3Int res, int zBegin, int zEnd);
+               static void setZeroBorder(float* field, Vec3Int res, int zBegin, int zEnd) {
+                       setZeroX(field, res, zBegin, zEnd);
+                       setZeroY(field, res, zBegin, zEnd);
+                       setZeroZ(field, res, zBegin, zEnd);
                };
 
+               
+
                // static advection functions, also used by WTURBULENCE
                static void advectFieldSemiLagrange(const float dt, const float* velx, const float* vely,  const float* velz,
-                               float* oldField, float* newField, Vec3Int res);
-               static void advectFieldMacCormack(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
-                               float* oldField, float* newField, float* temp1, float* temp2, Vec3Int res, const unsigned char* obstacles);
+                               float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd);
+               static void advectFieldMacCormack1(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
+                               float* oldField, float* tempResult, Vec3Int res, int zBegin, int zEnd);
+               static void advectFieldMacCormack2(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
+                               float* oldField, float* newField, float* tempResult, float* temp1,Vec3Int res, const unsigned char* obstacles, int zBegin, int zEnd);
+
+
+               // temp ones for testing
+               /*static void advectFieldMacCormack(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
+                               float* oldField, float* newField, float* temp1, float* temp2, Vec3Int res, const unsigned char* obstacles);*/
+               /*static void advectFieldSemiLagrange2(const float dt, const float* velx, const float* vely,  const float* velz,
+                               float* oldField, float* newField, Vec3Int res);*/
 
                // maccormack helper functions
                static void clampExtrema(const float dt, const float* xVelocity, const float* yVelocity,  const float* zVelocity,
-                               float* oldField, float* newField, Vec3Int res);
+                               float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd);
                static void clampOutsideRays(const float dt, const float* xVelocity, const float* yVelocity,  const float* zVelocity,
-                               float* oldField, float* newField, Vec3Int res, const unsigned char* obstacles, const float *oldAdvection);
+                               float* oldField, float* newField, Vec3Int res, const unsigned char* obstacles, const float *oldAdvection, int zBegin, int zEnd);
+
+
 
                // output helper functions
                // static void writeImageSliceXY(const float *field, Vec3Int res, int slice, string prefix, int picCnt, float scale=1.);
index 7d078d86d61de9be98213ca00885515d01f3327a..8b961494ce59456408e4dbe9bc9d0876826c2289 100644 (file)
 // FLUID_3D.cpp: implementation of the FLUID_3D class.
 //
 //////////////////////////////////////////////////////////////////////
+// Both solvers optimized by merging loops and precalculating
+// stuff used in iteration loop.
+//             - MiikaH
+//////////////////////////////////////////////////////////////////////
 
 #include "FLUID_3D.h"
 #include <cstring>
 #define SOLVER_ACCURACY 1e-06
 
-void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip)
+//////////////////////////////////////////////////////////////////////
+// solve the heat equation with CG
+//////////////////////////////////////////////////////////////////////
+void FLUID_3D::solveHeat(float* field, float* b, unsigned char* skip)
 {
        int x, y, z;
+       const int twoxr = 2 * _xRes;
        size_t index;
-       float *_q, *_Precond, *_h, *_residual, *_direction;
+       const float heatConst = _dt * _heatDiffusion / (_dx * _dx);
+       float *_q, *_residual, *_direction, *_Acenter;
 
        // i = 0
        int i = 0;
@@ -36,246 +45,201 @@ void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip)
        _residual     = new float[_totalCells]; // set 0
        _direction    = new float[_totalCells]; // set 0
        _q            = new float[_totalCells]; // set 0
-       _h                        = new float[_totalCells]; // set 0
-       _Precond          = new float[_totalCells]; // set 0
+       _Acenter       = new float[_totalCells]; // set 0
 
-       memset(_residual, 0, sizeof(float)*_xRes*_yRes*_zRes);
-       memset(_q, 0, sizeof(float)*_xRes*_yRes*_zRes);
-       memset(_direction, 0, sizeof(float)*_xRes*_yRes*_zRes);
-       memset(_h, 0, sizeof(float)*_xRes*_yRes*_zRes);
-       memset(_Precond, 0, sizeof(float)*_xRes*_yRes*_zRes);
-
-       // r = b - Ax
-       index = _slabSize + _xRes + 1;
-       for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-               for (y = 1; y < _yRes - 1; y++, index += 2)
-                 for (x = 1; x < _xRes - 1; x++, index++)
-                 {
-                       // if the cell is a variable
-                       float Acenter = 0.0f;
-                       if (!skip[index])
-                       {
-                         // set the matrix to the Poisson stencil in order
-                         if (!skip[index + 1]) Acenter += 1.;
-                         if (!skip[index - 1]) Acenter += 1.;
-                         if (!skip[index + _xRes]) Acenter += 1.;
-                         if (!skip[index - _xRes]) Acenter += 1.;
-                         if (!skip[index + _slabSize]) Acenter += 1.;
-                         if (!skip[index - _slabSize]) Acenter += 1.;
-                       }
-                   
-                       _residual[index] = b[index] - (Acenter * field[index] +  
-                         field[index - 1] * (skip[index - 1] ? 0.0 : -1.0f)+ 
-                         field[index + 1] * (skip[index + 1] ? 0.0 : -1.0f)+
-                         field[index - _xRes] * (skip[index - _xRes] ? 0.0 : -1.0f)+ 
-                         field[index + _xRes] * (skip[index + _xRes] ? 0.0 : -1.0f)+
-                         field[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -1.0f)+ 
-                         field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f) );
-                       _residual[index] = (skip[index]) ? 0.0f : _residual[index];
+       memset(_residual, 0, sizeof(float)*_totalCells);
+       memset(_q, 0, sizeof(float)*_totalCells);
+       memset(_direction, 0, sizeof(float)*_totalCells);
+       memset(_Acenter, 0, sizeof(float)*_totalCells);
 
-                       // P^-1
-                       if(Acenter < 1.0)
-                               _Precond[index] = 0.0;
-                       else
-                               _Precond[index] = 1.0 / Acenter;
+       float deltaNew = 0.0f;
 
-                       // p = P^-1 * r
-                       _direction[index] = _residual[index] * _Precond[index];
-                 }
+  // r = b - Ax
+  index = _slabSize + _xRes + 1;
+  for (z = 1; z < _zRes - 1; z++, index += twoxr)
+    for (y = 1; y < _yRes - 1; y++, index += 2)
+      for (x = 1; x < _xRes - 1; x++, index++)
+      {
+        // if the cell is a variable
+        _Acenter[index] = 1.0f;
+        if (!skip[index])
+        {
+          // set the matrix to the Poisson stencil in order
+          if (!skip[index + 1]) _Acenter[index] += heatConst;
+          if (!skip[index - 1]) _Acenter[index] += heatConst;
+          if (!skip[index + _xRes]) _Acenter[index] += heatConst;
+          if (!skip[index - _xRes]) _Acenter[index] += heatConst;
+          if (!skip[index + _slabSize]) _Acenter[index] += heatConst;
+          if (!skip[index - _slabSize]) _Acenter[index] += heatConst;
+
+                 _residual[index] = b[index] - (_Acenter[index] * field[index] + 
+          field[index - 1] * (skip[index - 1] ? 0.0 : -heatConst) + 
+          field[index + 1] * (skip[index + 1] ? 0.0 : -heatConst) +
+          field[index - _xRes] * (skip[index - _xRes] ? 0.0 : -heatConst) + 
+          field[index + _xRes] * (skip[index + _xRes] ? 0.0 : -heatConst) +
+          field[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -heatConst) + 
+          field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -heatConst));
+        }
+               else
+               {
+        _residual[index] = 0.0f;
+               }
 
-       // deltaNew = transpose(r) * p
-       float deltaNew = 0.0f;
-       index = _slabSize + _xRes + 1;
-       for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-               for (y = 1; y < _yRes - 1; y++, index += 2)
-                 for (x = 1; x < _xRes - 1; x++, index++)
-                       deltaNew += _residual[index] * _direction[index];
+               _direction[index] = _residual[index];
+               deltaNew += _residual[index] * _residual[index];
+      }
 
-       // delta0 = deltaNew
-       // float delta0 = deltaNew;
 
   // While deltaNew > (eps^2) * delta0
   const float eps  = SOLVER_ACCURACY;
-  //while ((i < _iterations) && (deltaNew > eps*delta0))
   float maxR = 2.0f * eps;
-  // while (i < _iterations)
-  while ((i < _iterations) && (maxR > 0.001*eps))
+  while ((i < _iterations) && (maxR > eps))
   {
-    // (s) q = Ad (p)
+    // q = Ad
+       float alpha = 0.0f;
+
     index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
+    for (z = 1; z < _zRes - 1; z++, index += twoxr)
       for (y = 1; y < _yRes - 1; y++, index += 2)
         for (x = 1; x < _xRes - 1; x++, index++)
         {
           // if the cell is a variable
-          float Acenter = 0.0f;
           if (!skip[index])
           {
-            // set the matrix to the Poisson stencil in order
-            if (!skip[index + 1]) Acenter += 1.;
-            if (!skip[index - 1]) Acenter += 1.;
-            if (!skip[index + _xRes]) Acenter += 1.;
-            if (!skip[index - _xRes]) Acenter += 1.;
-            if (!skip[index + _slabSize]) Acenter += 1.;
-            if (!skip[index - _slabSize]) Acenter += 1.;
+
+                       _q[index] = (_Acenter[index] * _direction[index] + 
+            _direction[index - 1] * (skip[index - 1] ? 0.0 : -heatConst) + 
+            _direction[index + 1] * (skip[index + 1] ? 0.0 : -heatConst) +
+            _direction[index - _xRes] * (skip[index - _xRes] ? 0.0 : -heatConst) + 
+            _direction[index + _xRes] * (skip[index + _xRes] ? 0.0 : -heatConst) +
+            _direction[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -heatConst) + 
+            _direction[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -heatConst));
           }
-          
-          _q[index] = Acenter * _direction[index] +  
-            _direction[index - 1] * (skip[index - 1] ? 0.0 : -1.0f) + 
-            _direction[index + 1] * (skip[index + 1] ? 0.0 : -1.0f) +
-            _direction[index - _xRes] * (skip[index - _xRes] ? 0.0 : -1.0f) + 
-            _direction[index + _xRes] * (skip[index + _xRes] ? 0.0 : -1.0f)+
-            _direction[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -1.0f) + 
-            _direction[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f);
-          _q[index] = (skip[index]) ? 0.0f : _q[index];
+                 else
+                 {
+          _q[index] = 0.0f;
+                 }
+                 alpha += _direction[index] * _q[index];
         }
 
-    // alpha = deltaNew / (transpose(d) * q)
-    float alpha = 0.0f;
-    index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-          alpha += _direction[index] * _q[index];
     if (fabs(alpha) > 0.0f)
       alpha = deltaNew / alpha;
+       
+       float deltaOld = deltaNew;
+       deltaNew = 0.0f;
 
-    // x = x + alpha * d
-    index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-          field[index] += alpha * _direction[index];
+       maxR = 0.0f;
 
-    // r = r - alpha * q
-       maxR = 0.0;
     index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-        {
-          _residual[index] -= alpha * _q[index];
-                 // maxR = (_residual[index] > maxR) ? _residual[index] : maxR;
-        }
-
-       // if(maxR <= eps)
-       //      break;
-
-       // h = P^-1 * r
-        index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
+    for (z = 1; z < _zRes - 1; z++, index += twoxr)
       for (y = 1; y < _yRes - 1; y++, index += 2)
         for (x = 1; x < _xRes - 1; x++, index++)
                {
-                       _h[index] = _Precond[index] * _residual[index];
-               }
+          field[index] += alpha * _direction[index];
 
-    // deltaOld = deltaNew
-    float deltaOld = deltaNew;
+                 _residual[index] -= alpha * _q[index];
+          maxR = (_residual[index] > maxR) ? _residual[index] : maxR;
 
-    // deltaNew = transpose(r) * h
-    deltaNew = 0.0f;
-    index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-               {
-          deltaNew += _residual[index] * _h[index];
-                 maxR = (_residual[index]* _h[index] > maxR) ? _residual[index]* _h[index] : maxR;
+                 deltaNew += _residual[index] * _residual[index];
                }
 
-    // beta = deltaNew / deltaOld
     float beta = deltaNew / deltaOld;
 
-    // d = h + beta * d
     index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
+    for (z = 1; z < _zRes - 1; z++, index += twoxr)
       for (y = 1; y < _yRes - 1; y++, index += 2)
         for (x = 1; x < _xRes - 1; x++, index++)
-          _direction[index] = _h[index] + beta * _direction[index];
+         _direction[index] = _residual[index] + beta * _direction[index];
 
-    // i = i + 1
+       
     i++;
   }
-  // cout << i << " iterations converged to " << sqrt(maxR) << endl;
+  // cout << i << " iterations converged to " << maxR << endl;
 
-       if (_h) delete[] _h;
-       if (_Precond) delete[] _Precond;
        if (_residual) delete[] _residual;
        if (_direction) delete[] _direction;
        if (_q)       delete[] _q;
+       if (_Acenter)  delete[] _Acenter;
 }
 
-//////////////////////////////////////////////////////////////////////
-// solve the poisson equation with CG
-//////////////////////////////////////////////////////////////////////
-#if 0
-void FLUID_3D::solvePressure(float* field, float* b, unsigned char* skip)
+
+void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip)
 {
        int x, y, z;
        size_t index;
+       float *_q, *_Precond, *_h, *_residual, *_direction;
 
        // i = 0
        int i = 0;
 
+       _residual     = new float[_totalCells]; // set 0
+       _direction    = new float[_totalCells]; // set 0
+       _q            = new float[_totalCells]; // set 0
+       _h                        = new float[_totalCells]; // set 0
+       _Precond          = new float[_totalCells]; // set 0
+
        memset(_residual, 0, sizeof(float)*_xRes*_yRes*_zRes);
        memset(_q, 0, sizeof(float)*_xRes*_yRes*_zRes);
        memset(_direction, 0, sizeof(float)*_xRes*_yRes*_zRes);
+       memset(_h, 0, sizeof(float)*_xRes*_yRes*_zRes);
+       memset(_Precond, 0, sizeof(float)*_xRes*_yRes*_zRes);
 
-  // r = b - Ax
-  index = _slabSize + _xRes + 1;
-  for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-    for (y = 1; y < _yRes - 1; y++, index += 2)
-      for (x = 1; x < _xRes - 1; x++, index++)
-      {
-        // if the cell is a variable
-        float Acenter = 0.0f;
-        if (!skip[index])
-        {
-          // set the matrix to the Poisson stencil in order
-          if (!skip[index + 1]) Acenter += 1.;
-          if (!skip[index - 1]) Acenter += 1.;
-          if (!skip[index + _xRes]) Acenter += 1.;
-          if (!skip[index - _xRes]) Acenter += 1.;
-          if (!skip[index + _slabSize]) Acenter += 1.;
-          if (!skip[index - _slabSize]) Acenter += 1.;
-        }
-        
-        _residual[index] = b[index] - (Acenter * field[index] +  
-          field[index - 1] * (skip[index - 1] ? 0.0 : -1.0f)+ 
-          field[index + 1] * (skip[index + 1] ? 0.0 : -1.0f)+
-          field[index - _xRes] * (skip[index - _xRes] ? 0.0 : -1.0f)+ 
-          field[index + _xRes] * (skip[index + _xRes] ? 0.0 : -1.0f)+
-          field[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -1.0f)+ 
-          field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f) );
-        _residual[index] = (skip[index]) ? 0.0f : _residual[index];
-      }
-       
+       float deltaNew = 0.0f;
 
-  // d = r
-  index = _slabSize + _xRes + 1;
-  for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-    for (y = 1; y < _yRes - 1; y++, index += 2)
-      for (x = 1; x < _xRes - 1; x++, index++)
-        _direction[index] = _residual[index];
+       // r = b - Ax
+       index = _slabSize + _xRes + 1;
+       for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
+               for (y = 1; y < _yRes - 1; y++, index += 2)
+                 for (x = 1; x < _xRes - 1; x++, index++)
+                 {
+                       // if the cell is a variable
+                       float Acenter = 0.0f;
+                       if (!skip[index])
+                       {
+                         // set the matrix to the Poisson stencil in order
+                         if (!skip[index + 1]) Acenter += 1.;
+                         if (!skip[index - 1]) Acenter += 1.;
+                         if (!skip[index + _xRes]) Acenter += 1.;
+                         if (!skip[index - _xRes]) Acenter += 1.;
+                         if (!skip[index + _slabSize]) Acenter += 1.;
+                         if (!skip[index - _slabSize]) Acenter += 1.;
 
-  // deltaNew = transpose(r) * r
-  float deltaNew = 0.0f;
-  index = _slabSize + _xRes + 1;
-  for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-    for (y = 1; y < _yRes - 1; y++, index += 2)
-      for (x = 1; x < _xRes - 1; x++, index++)
-        deltaNew += _residual[index] * _residual[index];
+                         _residual[index] = b[index] - (Acenter * field[index] +  
+                         field[index - 1] * (skip[index - 1] ? 0.0 : -1.0f)+ 
+                         field[index + 1] * (skip[index + 1] ? 0.0 : -1.0f)+
+                         field[index - _xRes] * (skip[index - _xRes] ? 0.0 : -1.0f)+ 
+                         field[index + _xRes] * (skip[index + _xRes] ? 0.0 : -1.0f)+
+                         field[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -1.0f)+ 
+                         field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f) );
+                       }
+                       else
+                       {
+                       _residual[index] = 0.0f;
+                       }
+
+                       // P^-1
+                       if(Acenter < 1.0)
+                               _Precond[index] = 0.0;
+                       else
+                               _Precond[index] = 1.0 / Acenter;
+
+                       // p = P^-1 * r
+                       _direction[index] = _residual[index] * _Precond[index];
+
+                       deltaNew += _residual[index] * _direction[index];
+                 }
 
-  // delta0 = deltaNew
-  // float delta0 = deltaNew;
 
   // While deltaNew > (eps^2) * delta0
   const float eps  = SOLVER_ACCURACY;
+  //while ((i < _iterations) && (deltaNew > eps*delta0))
   float maxR = 2.0f * eps;
-  while ((i < _iterations) && (maxR > eps))
+  // while (i < _iterations)
+  while ((i < _iterations) && (maxR > 0.001*eps))
   {
-    // q = Ad
+
+       float alpha = 0.0f;
+
     index = _slabSize + _xRes + 1;
     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
       for (y = 1; y < _yRes - 1; y++, index += 2)
@@ -292,233 +256,71 @@ void FLUID_3D::solvePressure(float* field, float* b, unsigned char* skip)
             if (!skip[index - _xRes]) Acenter += 1.;
             if (!skip[index + _slabSize]) Acenter += 1.;
             if (!skip[index - _slabSize]) Acenter += 1.;
-          }
-          
-          _q[index] = Acenter * _direction[index] +  
+
+                       _q[index] = Acenter * _direction[index] +  
             _direction[index - 1] * (skip[index - 1] ? 0.0 : -1.0f) + 
             _direction[index + 1] * (skip[index + 1] ? 0.0 : -1.0f) +
             _direction[index - _xRes] * (skip[index - _xRes] ? 0.0 : -1.0f) + 
             _direction[index + _xRes] * (skip[index + _xRes] ? 0.0 : -1.0f)+
             _direction[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -1.0f) + 
             _direction[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f);
-          _q[index] = (skip[index]) ? 0.0f : _q[index];
-        }
-
-    // alpha = deltaNew / (transpose(d) * q)
-    float alpha = 0.0f;
-    index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-          alpha += _direction[index] * _q[index];
-    if (fabs(alpha) > 0.0f)
-      alpha = deltaNew / alpha;
-
-    // x = x + alpha * d
-    index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-          field[index] += alpha * _direction[index];
-
-    // r = r - alpha * q
-    maxR = 0.0f;
-    index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-        {
-          _residual[index] -= alpha * _q[index];
-          maxR = (_residual[index] > maxR) ? _residual[index] : maxR;
-        }
-
-    // deltaOld = deltaNew
-    float deltaOld = deltaNew;
-
-    // deltaNew = transpose(r) * r
-    deltaNew = 0.0f;
-    index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-          deltaNew += _residual[index] * _residual[index];
-
-    // beta = deltaNew / deltaOld
-    float beta = deltaNew / deltaOld;
-
-    // d = r + beta * d
-    index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-          _direction[index] = _residual[index] + beta * _direction[index];
-
-    // i = i + 1
-    i++;
-  }
-  // cout << i << " iterations converged to " << maxR << endl;
-}
-#endif
-
-//////////////////////////////////////////////////////////////////////
-// solve the heat equation with CG
-//////////////////////////////////////////////////////////////////////
-void FLUID_3D::solveHeat(float* field, float* b, unsigned char* skip)
-{
-       int x, y, z;
-       size_t index;
-       const float heatConst = _dt * _heatDiffusion / (_dx * _dx);
-       float *_q, *_residual, *_direction;
-
-       // i = 0
-       int i = 0;
-
-       _residual     = new float[_totalCells]; // set 0
-       _direction    = new float[_totalCells]; // set 0
-       _q            = new float[_totalCells]; // set 0
-
-       memset(_residual, 0, sizeof(float)*_xRes*_yRes*_zRes);
-       memset(_q, 0, sizeof(float)*_xRes*_yRes*_zRes);
-       memset(_direction, 0, sizeof(float)*_xRes*_yRes*_zRes);
+          }
+                 else
+                 {
+          _q[index] = 0.0f;
+                 }
 
-  // r = b - Ax
-  index = _slabSize + _xRes + 1;
-  for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-    for (y = 1; y < _yRes - 1; y++, index += 2)
-      for (x = 1; x < _xRes - 1; x++, index++)
-      {
-        // if the cell is a variable
-        float Acenter = 1.0f;
-        if (!skip[index])
-        {
-          // set the matrix to the Poisson stencil in order
-          if (!skip[index + 1]) Acenter += heatConst;
-          if (!skip[index - 1]) Acenter += heatConst;
-          if (!skip[index + _xRes]) Acenter += heatConst;
-          if (!skip[index - _xRes]) Acenter += heatConst;
-          if (!skip[index + _slabSize]) Acenter += heatConst;
-          if (!skip[index - _slabSize]) Acenter += heatConst;
+                 alpha += _direction[index] * _q[index];
         }
-        
-        _residual[index] = b[index] - (Acenter * field[index] + 
-          field[index - 1] * (skip[index - 1] ? 0.0 : -heatConst) + 
-          field[index + 1] * (skip[index + 1] ? 0.0 : -heatConst) +
-          field[index - _xRes] * (skip[index - _xRes] ? 0.0 : -heatConst) + 
-          field[index + _xRes] * (skip[index + _xRes] ? 0.0 : -heatConst) +
-          field[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -heatConst) + 
-          field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -heatConst));
-        _residual[index] = (skip[index]) ? 0.0f : _residual[index];
-      }
 
-  // d = r
-  index = _slabSize + _xRes + 1;
-  for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-    for (y = 1; y < _yRes - 1; y++, index += 2)
-      for (x = 1; x < _xRes - 1; x++, index++)
-        _direction[index] = _residual[index];
 
-  // deltaNew = transpose(r) * r
-  float deltaNew = 0.0f;
-  index = _slabSize + _xRes + 1;
-  for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-    for (y = 1; y < _yRes - 1; y++, index += 2)
-      for (x = 1; x < _xRes - 1; x++, index++)
-        deltaNew += _residual[index] * _residual[index];
+    if (fabs(alpha) > 0.0f)
+      alpha = deltaNew / alpha;
 
-  // delta0 = deltaNew
-  // float delta0 = deltaNew;
+       float deltaOld = deltaNew;
+       deltaNew = 0.0f;
 
-  // While deltaNew > (eps^2) * delta0
-  const float eps  = SOLVER_ACCURACY;
-  float maxR = 2.0f * eps;
-  while ((i < _iterations) && (maxR > eps))
-  {
-    // q = Ad
-    index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-        {
-          // if the cell is a variable
-          float Acenter = 1.0f;
-          if (!skip[index])
-          {
-            // set the matrix to the Poisson stencil in order
-            if (!skip[index + 1]) Acenter += heatConst;
-            if (!skip[index - 1]) Acenter += heatConst;
-            if (!skip[index + _xRes]) Acenter += heatConst;
-            if (!skip[index - _xRes]) Acenter += heatConst;
-            if (!skip[index + _slabSize]) Acenter += heatConst;
-            if (!skip[index - _slabSize]) Acenter += heatConst;
-          }
-
-          _q[index] = (Acenter * _direction[index] + 
-            _direction[index - 1] * (skip[index - 1] ? 0.0 : -heatConst) + 
-            _direction[index + 1] * (skip[index + 1] ? 0.0 : -heatConst) +
-            _direction[index - _xRes] * (skip[index - _xRes] ? 0.0 : -heatConst) + 
-            _direction[index + _xRes] * (skip[index + _xRes] ? 0.0 : -heatConst) +
-            _direction[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -heatConst) + 
-            _direction[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -heatConst));
-         
-          _q[index] = (skip[index]) ? 0.0f : _q[index];
-        }
+       maxR = 0.0;
 
-    // alpha = deltaNew / (transpose(d) * q)
-    float alpha = 0.0f;
-    index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-          alpha += _direction[index] * _q[index];
-    if (fabs(alpha) > 0.0f)
-      alpha = deltaNew / alpha;
+       float tmp;
 
     // x = x + alpha * d
     index = _slabSize + _xRes + 1;
     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
       for (y = 1; y < _yRes - 1; y++, index += 2)
         for (x = 1; x < _xRes - 1; x++, index++)
+               {
           field[index] += alpha * _direction[index];
 
-    // r = r - alpha * q
-    maxR = 0.0f;
-    index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-        {
-          _residual[index] -= alpha * _q[index];
-          maxR = (_residual[index] > maxR) ? _residual[index] : maxR;
-        }
+                 _residual[index] -= alpha * _q[index];
 
-    // deltaOld = deltaNew
-    float deltaOld = deltaNew;
+                 _h[index] = _Precond[index] * _residual[index];
+
+                 tmp = _residual[index] * _h[index];
+                 deltaNew += tmp;
+                 maxR = (tmp > maxR) ? tmp : maxR;
+
+               }
 
-    // deltaNew = transpose(r) * r
-    deltaNew = 0.0f;
-    index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-          deltaNew += _residual[index] * _residual[index];
 
     // beta = deltaNew / deltaOld
     float beta = deltaNew / deltaOld;
 
-    // d = r + beta * d
+    // d = h + beta * d
     index = _slabSize + _xRes + 1;
     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
       for (y = 1; y < _yRes - 1; y++, index += 2)
         for (x = 1; x < _xRes - 1; x++, index++)
-          _direction[index] = _residual[index] + beta * _direction[index];
+          _direction[index] = _h[index] + beta * _direction[index];
 
     // i = i + 1
     i++;
   }
-  // cout << i << " iterations converged to " << maxR << endl;
+  // cout << i << " iterations converged to " << sqrt(maxR) << endl;
 
+       if (_h) delete[] _h;
+       if (_Precond) delete[] _Precond;
        if (_residual) delete[] _residual;
        if (_direction) delete[] _direction;
        if (_q)       delete[] _q;
 }
-
index afeca2b1faaea6438d4fcfc3eb7dad9e4f0980b4..89ac2c74e15185992cc6c22c43a4d23e64f5c80b 100644 (file)
 // FLUID_3D.cpp: implementation of the static functions of the FLUID_3D class.
 //
 //////////////////////////////////////////////////////////////////////
+// Heavy parallel optimization done. Many of the old functions now
+// take begin and end parameters and process only specified part of the data.
+// Some functions were divided into multiple ones.
+//             - MiikaH
+//////////////////////////////////////////////////////////////////////
 
 #include <zlib.h>
 #include "FLUID_3D.h"
@@ -75,11 +80,11 @@ void FLUID_3D::addSmokeTestCase(float* field, Vec3Int res)
 //////////////////////////////////////////////////////////////////////
 // set x direction to Neumann boundary conditions
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::setNeumannX(float* field, Vec3Int res)
+void FLUID_3D::setNeumannX(float* field, Vec3Int res, int zBegin, int zEnd)
 {
        const int slabSize = res[0] * res[1];
        int index;
-       for (int z = 0; z < res[2]; z++)
+       for (int z = zBegin; z < zEnd; z++)
                for (int y = 0; y < res[1]; y++)
                {
                        // left slab
@@ -93,7 +98,7 @@ void FLUID_3D::setNeumannX(float* field, Vec3Int res)
 
        // fix, force top slab to only allow outwards flux
        for (int y = 0; y < res[1]; y++)
-               for (int z = 0; z < res[2]; z++)
+               for (int z = zBegin; z < zEnd; z++)
                {
                        // top slab
                        index = y * res[0] + z * slabSize;
@@ -107,11 +112,11 @@ void FLUID_3D::setNeumannX(float* field, Vec3Int res)
 //////////////////////////////////////////////////////////////////////
 // set y direction to Neumann boundary conditions
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::setNeumannY(float* field, Vec3Int res)
+void FLUID_3D::setNeumannY(float* field, Vec3Int res, int zBegin, int zEnd)
 {
        const int slabSize = res[0] * res[1];
        int index;
-       for (int z = 0; z < res[2]; z++)
+       for (int z = zBegin; z < zEnd; z++)
                for (int x = 0; x < res[0]; x++)
                {
                        // bottom slab
@@ -124,7 +129,7 @@ void FLUID_3D::setNeumannY(float* field, Vec3Int res)
                }
 
        // fix, force top slab to only allow outwards flux
-       for (int z = 0; z < res[2]; z++)
+       for (int z = zBegin; z < zEnd; z++)
                for (int x = 0; x < res[0]; x++)
                {
                        // top slab
@@ -140,22 +145,36 @@ void FLUID_3D::setNeumannY(float* field, Vec3Int res)
 //////////////////////////////////////////////////////////////////////
 // set z direction to Neumann boundary conditions
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::setNeumannZ(float* field, Vec3Int res)
+void FLUID_3D::setNeumannZ(float* field, Vec3Int res, int zBegin, int zEnd)
 {
        const int slabSize = res[0] * res[1];
        const int totalCells = res[0] * res[1] * res[2];
+       const int cellsslab = totalCells - slabSize;
        int index;
+
+       index = 0;
+       if (zBegin == 0)
        for (int y = 0; y < res[1]; y++)
-               for (int x = 0; x < res[0]; x++)
+               for (int x = 0; x < res[0]; x++, index++)
                {
                        // front slab
-                       index = x + y * res[0];
                        field[index] = field[index + 2 * slabSize];
+               }
+
+       if (zEnd == res[2])
+       {
+       index = 0;
+       int indexx = 0;
+
+       for (int y = 0; y < res[1]; y++)
+               for (int x = 0; x < res[0]; x++, index++)
+               {
 
                        // back slab
-                       index += totalCells - slabSize;
-                       field[index] = field[index - 2 * slabSize];
+                       indexx = index + cellsslab;
+                       field[indexx] = field[indexx - 2 * slabSize];
                }
+       
 
        // fix, force top slab to only allow outwards flux
        for (int y = 0; y < res[1]; y++)
@@ -163,22 +182,24 @@ void FLUID_3D::setNeumannZ(float* field, Vec3Int res)
                {
                        // top slab
                        index = x + y * res[0];
-                       index += totalCells - slabSize;
+                       index += cellsslab;
                        if(field[index]<0.) field[index] = 0.;
                        index -= slabSize;
                        if(field[index]<0.) field[index] = 0.;
                }
+
+       }       // zEnd == res[2]
                
 }
 
 //////////////////////////////////////////////////////////////////////
 // set x direction to zero
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::setZeroX(float* field, Vec3Int res)
+void FLUID_3D::setZeroX(float* field, Vec3Int res, int zBegin, int zEnd)
 {
        const int slabSize = res[0] * res[1];
        int index;
-       for (int z = 0; z < res[2]; z++)
+       for (int z = zBegin; z < zEnd; z++)
                for (int y = 0; y < res[1]; y++)
                {
                        // left slab
@@ -194,11 +215,11 @@ void FLUID_3D::setZeroX(float* field, Vec3Int res)
 //////////////////////////////////////////////////////////////////////
 // set y direction to zero
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::setZeroY(float* field, Vec3Int res)
+void FLUID_3D::setZeroY(float* field, Vec3Int res, int zBegin, int zEnd)
 {
        const int slabSize = res[0] * res[1];
        int index;
-       for (int z = 0; z < res[2]; z++)
+       for (int z = zBegin; z < zEnd; z++)
                for (int x = 0; x < res[0]; x++)
                {
                        // bottom slab
@@ -214,32 +235,44 @@ void FLUID_3D::setZeroY(float* field, Vec3Int res)
 //////////////////////////////////////////////////////////////////////
 // set z direction to zero
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::setZeroZ(float* field, Vec3Int res)
+void FLUID_3D::setZeroZ(float* field, Vec3Int res, int zBegin, int zEnd)
 {
        const int slabSize = res[0] * res[1];
        const int totalCells = res[0] * res[1] * res[2];
-       int index;
+
+       int index = 0;
+       if ((zBegin == 0))
        for (int y = 0; y < res[1]; y++)
-               for (int x = 0; x < res[0]; x++)
+               for (int x = 0; x < res[0]; x++, index++)
                {
                        // front slab
-                       index = x + y * res[0];
                        field[index] = 0.0f;
+    }
 
-                       // back slab
-                       index += totalCells - slabSize;
-                       field[index] = 0.0f;
-               }
- }
+       if (zEnd == res[2])
+       {
+               index=0;
+               int indexx=0;
+               const int cellsslab = totalCells - slabSize;
+
+               for (int y = 0; y < res[1]; y++)
+                       for (int x = 0; x < res[0]; x++, index++)
+                       {
 
+                               // back slab
+                               indexx = index + cellsslab;
+                               field[indexx] = 0.0f;
+                       }
+       }
+ }
 //////////////////////////////////////////////////////////////////////
 // copy grid boundary
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::copyBorderX(float* field, Vec3Int res)
+void FLUID_3D::copyBorderX(float* field, Vec3Int res, int zBegin, int zEnd)
 {
        const int slabSize = res[0] * res[1];
        int index;
-       for (int z = 0; z < res[2]; z++)
+       for (int z = zBegin; z < zEnd; z++)
                for (int y = 0; y < res[1]; y++)
                {
                        // left slab
@@ -251,12 +284,12 @@ void FLUID_3D::copyBorderX(float* field, Vec3Int res)
                        field[index] = field[index - 1];
                }
 }
-void FLUID_3D::copyBorderY(float* field, Vec3Int res)
+void FLUID_3D::copyBorderY(float* field, Vec3Int res, int zBegin, int zEnd)
 {
        const int slabSize = res[0] * res[1];
-       const int totalCells = res[0] * res[1] * res[2];
+       //const int totalCells = res[0] * res[1] * res[2];
        int index;
-       for (int z = 0; z < res[2]; z++)
+       for (int z = zBegin; z < zEnd; z++)
                for (int x = 0; x < res[0]; x++)
                {
                        // bottom slab
@@ -267,40 +300,49 @@ void FLUID_3D::copyBorderY(float* field, Vec3Int res)
                        field[index] = field[index - res[0]];
                }
 }
-void FLUID_3D::copyBorderZ(float* field, Vec3Int res)
+void FLUID_3D::copyBorderZ(float* field, Vec3Int res, int zBegin, int zEnd)
 {
        const int slabSize = res[0] * res[1];
        const int totalCells = res[0] * res[1] * res[2];
-       int index;
+       int index=0;
+
+       if ((zBegin == 0))
        for (int y = 0; y < res[1]; y++)
-               for (int x = 0; x < res[0]; x++)
+               for (int x = 0; x < res[0]; x++, index++)
                {
-                       // front slab
-                       index = x + y * res[0];
                        field[index] = field[index + slabSize]; 
+               }
+
+       if ((zEnd == res[2]))
+       {
+
+       index=0;
+       int indexx=0;
+       const int cellsslab = totalCells - slabSize;
+
+       for (int y = 0; y < res[1]; y++)
+               for (int x = 0; x < res[0]; x++, index++)
+               {
                        // back slab
-                       index += totalCells - slabSize;
-                       field[index] = field[index - slabSize];
+                       indexx = index + cellsslab;
+                       field[indexx] = field[indexx - slabSize];
                }
+       }
 }
 
 /////////////////////////////////////////////////////////////////////
 // advect field with the semi lagrangian method
 //////////////////////////////////////////////////////////////////////
 void FLUID_3D::advectFieldSemiLagrange(const float dt, const float* velx, const float* vely,  const float* velz,
-               float* oldField, float* newField, Vec3Int res)
+               float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd)
 {
        const int xres = res[0];
        const int yres = res[1];
        const int zres = res[2];
        const int slabSize = res[0] * res[1];
 
-       // scale dt up to grid resolution
-#if PARALLEL==1 && !_WIN32
-#pragma omp parallel
-#pragma omp for  schedule(static)
-#endif
-       for (int z = 0; z < zres; z++)
+
+       for (int z = zBegin; z < zEnd; z++)
                for (int y = 0; y < yres; y++)
                        for (int x = 0; x < xres; x++)
                        {
@@ -357,50 +399,69 @@ void FLUID_3D::advectFieldSemiLagrange(const float dt, const float* velx, const
                        }
 }
 
+
 /////////////////////////////////////////////////////////////////////
 // advect field with the maccormack method
 //
 // comments are the pseudocode from selle's paper
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::advectFieldMacCormack(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
-                               float* oldField, float* newField, float* temp1, float* temp2, Vec3Int res, const unsigned char* obstacles)
+void FLUID_3D::advectFieldMacCormack1(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
+                               float* oldField, float* tempResult, Vec3Int res, int zBegin, int zEnd)
 {
-       float* phiHatN  = temp1;
-       float* phiHatN1 = temp2;
        const int sx= res[0];
        const int sy= res[1];
        const int sz= res[2];
 
-       for (int x = 0; x < sx * sy * sz; x++)
-               phiHatN[x] = phiHatN1[x] = oldField[x];
+       /*for (int x = 0; x < sx * sy * sz; x++)
+               phiHatN[x] = phiHatN1[x] = oldField[x];*/       // not needed as all the values are written first
 
        float*& phiN    = oldField;
-       float*& phiN1   = newField;
+       float*& phiN1   = tempResult;
+
+
 
        // phiHatN1 = A(phiN)
-       advectFieldSemiLagrange(  dt, xVelocity, yVelocity, zVelocity, phiN, phiHatN1, res);
+       advectFieldSemiLagrange(  dt, xVelocity, yVelocity, zVelocity, phiN, phiN1, res, zBegin, zEnd);         // uses wide data from old field and velocities (both are whole)
+}
+
+
+
+void FLUID_3D::advectFieldMacCormack2(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
+                               float* oldField, float* newField, float* tempResult, float* temp1, Vec3Int res, const unsigned char* obstacles, int zBegin, int zEnd)
+{
+       float* phiHatN  = tempResult;
+       float* t1  = temp1;
+       const int sx= res[0];
+       const int sy= res[1];
+       const int sz= res[2];
+
+       float*& phiN    = oldField;
+       float*& phiN1   = newField;
+
+
 
        // phiHatN = A^R(phiHatN1)
-       advectFieldSemiLagrange( -1.0*dt, xVelocity, yVelocity, zVelocity, phiHatN1, phiHatN, res);
+       advectFieldSemiLagrange( -1.0*dt, xVelocity, yVelocity, zVelocity, phiHatN, t1, res, zBegin, zEnd);             // uses wide data from old field and velocities (both are whole)
 
        // phiN1 = phiHatN1 + (phiN - phiHatN) / 2
        const int border = 0; 
-       for (int z = border; z < sz-border; z++)
+       for (int z = zBegin+border; z < zEnd-border; z++)
                for (int y = border; y < sy-border; y++)
                        for (int x = border; x < sx-border; x++) {
                                int index = x + y * sx + z * sx*sy;
-                               phiN1[index] = phiHatN1[index] + (phiN[index] - phiHatN[index]) * 0.50f;
+                               phiN1[index] = phiHatN[index] + (phiN[index] - t1[index]) * 0.50f;
                                //phiN1[index] = phiHatN1[index]; // debug, correction off
                        }
-       copyBorderX(phiN1, res);
-       copyBorderY(phiN1, res);
-       copyBorderZ(phiN1, res);
+       copyBorderX(phiN1, res, zBegin, zEnd);
+       copyBorderY(phiN1, res, zBegin, zEnd);
+       copyBorderZ(phiN1, res, zBegin, zEnd);
 
        // clamp any newly created extrema
-       clampExtrema(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res);
+       clampExtrema(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, zBegin, zEnd);               // uses wide data from old field and velocities (both are whole)
 
        // if the error estimate was bad, revert to first order
-       clampOutsideRays(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, obstacles, phiHatN1);
+       clampOutsideRays(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, obstacles, phiHatN, zBegin, zEnd);       // phiHatN is only used at cells within thread range, so its ok
+
 } 
 
 
@@ -408,13 +469,21 @@ void FLUID_3D::advectFieldMacCormack(const float dt, const float* xVelocity, con
 // Clamp the extrema generated by the BFECC error correction
 //////////////////////////////////////////////////////////////////////
 void FLUID_3D::clampExtrema(const float dt, const float* velx, const float* vely,  const float* velz,
-               float* oldField, float* newField, Vec3Int res)
+               float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd)
 {
        const int xres= res[0];
        const int yres= res[1];
        const int zres= res[2];
        const int slabSize = res[0] * res[1];
-       for (int z = 1; z < zres-1; z++)
+
+       int bb=0;
+       int bt=0;
+
+       if (zBegin == 0) {bb = 1;}
+       if (zEnd == res[2]) {bt = 1;}
+
+
+       for (int z = zBegin+bb; z < zEnd-bt; z++)
                for (int y = 1; y < yres-1; y++)
                        for (int x = 1; x < xres-1; x++)
                        {
@@ -484,14 +553,20 @@ void FLUID_3D::clampExtrema(const float dt, const float* velx, const float* vely
 // incorrect
 //////////////////////////////////////////////////////////////////////
 void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* vely,  const float* velz,
-                               float* oldField, float* newField, Vec3Int res, const unsigned char* obstacles, const float *oldAdvection)
+                               float* oldField, float* newField, Vec3Int res, const unsigned char* obstacles, const float *oldAdvection, int zBegin, int zEnd)
 {
        const int sx= res[0];
        const int sy= res[1];
        const int sz= res[2];
        const int slabSize = res[0] * res[1];
 
-       for (int z = 1; z < sz-1; z++)
+       int bb=0;
+       int bt=0;
+
+       if (zBegin == 0) {bb = 1;}
+       if (zEnd == res[2]) {bt = 1;}
+
+       for (int z = zBegin+bb; z < zEnd-bt; z++)
                for (int y = 1; y < sy-1; y++)
                        for (int x = 1; x < sx-1; x++)
                        {
@@ -607,4 +682,3 @@ void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float*
                                }
                        } // xyz
 }
-
diff --git a/intern/smoke/intern/LU_HELPER.cpp b/intern/smoke/intern/LU_HELPER.cpp
new file mode 100644 (file)
index 0000000..43357e6
--- /dev/null
@@ -0,0 +1,136 @@
+
+#include "LU_HELPER.h"
+
+int isNonsingular (sLU LU_) {
+      for (int j = 0; j < 3; j++) {
+         if (LU_.values[j][j] == 0)
+            return 0;
+      }
+      return 1;
+}
+
+sLU computeLU( float a[3][3])
+{
+       sLU result;
+       int m=3;
+       int n=3;
+
+       //float LU_[3][3];
+       for (int i = 0; i < m; i++) {
+                       result.piv[i] = i;
+                       for (int j = 0; j < n; j++) result.values[i][j]=a[i][j];
+      }
+
+      result.pivsign = 1;
+      //Real *LUrowi = 0;;
+      //Array1D<Real> LUcolj(m);
+         //float *LUrowi = 0;
+         float LUcolj[3];
+
+      // Outer loop.
+
+      for (int j = 0; j < n; j++) {
+
+         // Make a copy of the j-th column to localize references.
+
+         for (int i = 0; i < m; i++) {
+            LUcolj[i] = result.values[i][j];
+         }
+
+         // Apply previous transformations.
+
+         for (int i = 0; i < m; i++) {
+                        //float LUrowi[3];
+                        //LUrowi = result.values[i];
+
+            // Most of the time is spent in the following dot product.
+
+            int kmax = min(i,j);
+            double s = 0.0;
+            for (int k = 0; k < kmax; k++) {
+               s += result.values[i][k]*LUcolj[k];
+            }
+
+            result.values[i][j] = LUcolj[i] -= s;
+         }
+   
+         // Find pivot and exchange if necessary.
+
+         int p = j;
+         for (int i = j+1; i < m; i++) {
+            if (abs(LUcolj[i]) > abs(LUcolj[p])) {
+               p = i;
+            }
+         }
+         if (p != j) {
+                   int k=0;
+            for (k = 0; k < n; k++) {
+               double t = result.values[p][k]; 
+                          result.values[p][k] = result.values[j][k]; 
+                          result.values[j][k] = t;
+            }
+            k = result.piv[p]; 
+                       result.piv[p] = result.piv[j]; 
+                       result.piv[j] = k;
+            result.pivsign = -result.pivsign;
+         }
+
+         // Compute multipliers.
+         
+         if ((j < m) && (result.values[j][j] != 0.0)) {
+            for (int i = j+1; i < m; i++) {
+               result.values[i][j] /= result.values[j][j];
+            }
+         }
+      }
+
+         return result;
+}
+
+void solveLU3x3(sLU& A, float x[3], float b[3])
+{
+  //TNT::Array1D<float> jamaB = TNT::Array1D<float>(3, &b[0]);
+  //TNT::Array1D<float> jamaX = A.solve(jamaB);
+
+       
+  // Solve A, B
+
+       {
+      if (!isNonsingular(A)) {
+        x[0]=0.0f;
+               x[1]=0.0f;
+               x[2]=0.0f;
+               return;
+      }
+
+
+         //Array1D<Real> Ax = permute_copy(b, piv);
+         float Ax[3];
+
+    // permute copy: b , A.piv
+       {
+         for (int i = 0; i < 3; i++) 
+               Ax[i] = b[A.piv[i]];
+       }
+
+      // Solve L*Y = B(piv)
+      for (int k = 0; k < 3; k++) {
+         for (int i = k+1; i < 3; i++) {
+               Ax[i] -= Ax[k]*A.values[i][k];
+            }
+         }
+      
+         // Solve U*X = Y;
+      for (int k = 2; k >= 0; k--) {
+            Ax[k] /= A.values[k][k];
+               for (int i = 0; i < k; i++) 
+               Ax[i] -= Ax[k]*A.values[i][k];
+      }
+     
+
+               x[0] = Ax[0];
+               x[1] = Ax[1];
+               x[2] = Ax[2];
+      return;
+       }
+}
index b3f3c5a1cb4e8378f093ea3ec404f3f5a5a7ee98..e79b4ffa01b92d27c4dc6d7774f5b1a4c73e54a9 100644 (file)
 // Copyright 2008 Theodore Kim and Nils Thuerey
 // 
 //////////////////////////////////////////////////////////////////////
+// Modified to not require TNT matrix library anymore. It was very slow
+// when being run in parallel. Required TNT JAMA:LU libraries were
+// converted into independent functions.
+//             - MiikaH
+//////////////////////////////////////////////////////////////////////
 
 #ifndef LU_HELPER_H
 #define LU_HELPER_H
 
-//////////////////////////////////////////////////////////////////////
-// Helper function, compute eigenvalues of 3x3 matrix
-//////////////////////////////////////////////////////////////////////
+#include <cmath>
+#include <algorithm>
 
-#include "tnt/jama_lu.h"
+using namespace std;
 
 //////////////////////////////////////////////////////////////////////
-// LU decomposition of 3x3 non-symmetric matrix
+// Helper function, compute eigenvalues of 3x3 matrix
 //////////////////////////////////////////////////////////////////////
-JAMA::LU<float> inline computeLU3x3(
-               float a[3][3])
-{
-               TNT::Array2D<float> A = TNT::Array2D<float>(3,3, &a[0][0]);
-               JAMA::LU<float> jLU= JAMA::LU<float>(A);
-               return jLU;
-}
 
-//////////////////////////////////////////////////////////////////////
-// LU decomposition of 3x3 non-symmetric matrix
-//////////////////////////////////////////////////////////////////////
-void inline solveLU3x3(
-    JAMA::LU<float>& A,
-    float x[3],
-    float b[3])
+struct sLU
 {
-  TNT::Array1D<float> jamaB = TNT::Array1D<float>(3, &b[0]);
-  TNT::Array1D<float> jamaX = A.solve(jamaB);
+       float values[3][3];
+       int pivsign;
+       int piv[3];
+};
+
+
+int isNonsingular (sLU LU_);
+sLU computeLU( float a[3][3]);
+void solveLU3x3(sLU& A, float x[3], float b[3]);
+
 
-  x[0] = jamaX[0];
-  x[1] = jamaX[1];
-  x[2] = jamaX[2];
-}
 #endif
index 7ea4bde3884fde9e242cf6aa58fc9ad3d11c0e28..6d43bc95471b5a50e7e70f5043e5cefe5e5c7bcc 100644 (file)
 // 
 // WTURBULENCE handling
 ///////////////////////////////////////////////////////////////////////////////////
+// Parallelized turbulence even further. TNT matrix library functions
+// rewritten to improve performance.
+//             - MiikaH
+//////////////////////////////////////////////////////////////////////
 
 #include "WTURBULENCE.h"
 #include "INTERPOLATE.h"
@@ -29,6 +33,7 @@
 #include "LU_HELPER.h"
 #include "SPHERE.h"
 #include <zlib.h>
+#include <math.h>
 
 // needed to access static advection functions
 #include "FLUID_3D.h"
@@ -278,27 +283,34 @@ static float minDz(int x, int y, int z, float* input, Vec3Int res)
 // Beware -- uses big density maccormack as temporary arrays
 ////////////////////////////////////////////////////////////////////// 
 void WTURBULENCE::advectTextureCoordinates (float dtOrg, float* xvel, float* yvel, float* zvel, float *tempBig1, float *tempBig2) {
+
   // advection
   SWAP_POINTERS(_tcTemp, _tcU);
-  FLUID_3D::copyBorderX(_tcTemp, _resSm);
-  FLUID_3D::copyBorderY(_tcTemp, _resSm);
-  FLUID_3D::copyBorderZ(_tcTemp, _resSm);
-  FLUID_3D::advectFieldMacCormack(dtOrg, xvel, yvel, zvel, 
-      _tcTemp, _tcU, tempBig1, tempBig2, _resSm, NULL);
+  FLUID_3D::copyBorderX(_tcTemp, _resSm, 0 , _resSm[2]);
+  FLUID_3D::copyBorderY(_tcTemp, _resSm, 0 , _resSm[2]);
+  FLUID_3D::copyBorderZ(_tcTemp, _resSm, 0 , _resSm[2]);
+  FLUID_3D::advectFieldMacCormack1(dtOrg, xvel, yvel, zvel, 
+      _tcTemp, tempBig1, _resSm, 0 , _resSm[2]);
+  FLUID_3D::advectFieldMacCormack2(dtOrg, xvel, yvel, zvel, 
+      _tcTemp, _tcU, tempBig1, tempBig2, _resSm, NULL, 0 , _resSm[2]);
 
   SWAP_POINTERS(_tcTemp, _tcV);
-  FLUID_3D::copyBorderX(_tcTemp, _resSm);
-  FLUID_3D::copyBorderY(_tcTemp, _resSm);
-  FLUID_3D::copyBorderZ(_tcTemp, _resSm);
-  FLUID_3D::advectFieldMacCormack(dtOrg, xvel, yvel, zvel, 
-      _tcTemp, _tcV, tempBig1, tempBig2, _resSm, NULL);
+  FLUID_3D::copyBorderX(_tcTemp, _resSm, 0 , _resSm[2]);
+  FLUID_3D::copyBorderY(_tcTemp, _resSm, 0 , _resSm[2]);
+  FLUID_3D::copyBorderZ(_tcTemp, _resSm, 0 , _resSm[2]);
+  FLUID_3D::advectFieldMacCormack1(dtOrg, xvel, yvel, zvel, 
+      _tcTemp, tempBig1, _resSm, 0 , _resSm[2]);
+  FLUID_3D::advectFieldMacCormack2(dtOrg, xvel, yvel, zvel, 
+      _tcTemp, _tcV, tempBig1, tempBig2, _resSm, NULL, 0 , _resSm[2]);
 
   SWAP_POINTERS(_tcTemp, _tcW);
-  FLUID_3D::copyBorderX(_tcTemp, _resSm);
-  FLUID_3D::copyBorderY(_tcTemp, _resSm);
-  FLUID_3D::copyBorderZ(_tcTemp, _resSm);
-  FLUID_3D::advectFieldMacCormack(dtOrg, xvel, yvel, zvel, 
-      _tcTemp, _tcW, tempBig1, tempBig2, _resSm, NULL);
+  FLUID_3D::copyBorderX(_tcTemp, _resSm, 0 , _resSm[2]);
+  FLUID_3D::copyBorderY(_tcTemp, _resSm, 0 , _resSm[2]);
+  FLUID_3D::copyBorderZ(_tcTemp, _resSm, 0 , _resSm[2]);
+  FLUID_3D::advectFieldMacCormack1(dtOrg, xvel, yvel, zvel, 
+      _tcTemp, tempBig1, _resSm, 0 , _resSm[2]);
+  FLUID_3D::advectFieldMacCormack2(dtOrg, xvel, yvel, zvel, 
+      _tcTemp, _tcW, tempBig1, tempBig2, _resSm, NULL, 0 , _resSm[2]);
 }
 
 //////////////////////////////////////////////////////////////////////
@@ -325,9 +337,9 @@ void WTURBULENCE::computeEigenvalues(float *_eigMin, float *_eigMax) {
 
         // ONLY compute the eigenvalues after checking that the matrix
         // is nonsingular
-        JAMA::LU<float> LU = computeLU3x3(jacobian);
+        sLU LU = computeLU(jacobian);
 
-        if (LU.isNonsingular())
+        if (isNonsingular(LU))
         {
           // get the analytic eigenvalues, quite slow right now...
           Vec3 eigenvalues = Vec3(1.);
@@ -422,9 +434,9 @@ void WTURBULENCE::computeEnergy(float *_energy, float* xvel, float* yvel, float*
   for (int x = 0; x < _totalCellsSm; x++) 
     _energy[x] = 0.5f * (xvel[x] * xvel[x] + yvel[x] * yvel[x] + zvel[x] * zvel[x]);
 
-  FLUID_3D::copyBorderX(_energy, _resSm);
-  FLUID_3D::copyBorderY(_energy, _resSm);
-  FLUID_3D::copyBorderZ(_energy, _resSm);
+  FLUID_3D::copyBorderX(_energy, _resSm, 0 , _resSm[2]);
+  FLUID_3D::copyBorderY(_energy, _resSm, 0 , _resSm[2]);
+  FLUID_3D::copyBorderZ(_energy, _resSm, 0 , _resSm[2]);
 
   // pseudo-march the values into the obstacles
   // the wavelet upsampler only uses a 3x3 support neighborhood, so
@@ -563,7 +575,7 @@ Vec3 WTURBULENCE::WVelocityWithJacobian(Vec3 orgPos, float* xUnwarped, float* yU
 //////////////////////////////////////////////////////////////////////
 // perform an actual noise advection step
 //////////////////////////////////////////////////////////////////////
-void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel, float* zvel, unsigned char *obstacles) 
+/*void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel, float* zvel, unsigned char *obstacles) 
 {
        // enlarge timestep to match grid
        const float dt = dtOrg * _amplify;
@@ -689,9 +701,9 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel,
   const float dtSubdiv = dt / (float)totalSubsteps;
 
   // set boundaries of big velocity grid
-  FLUID_3D::setZeroX(bigUx, _resBig); 
-  FLUID_3D::setZeroY(bigUy, _resBig); 
-  FLUID_3D::setZeroZ(bigUz, _resBig);
+  FLUID_3D::setZeroX(bigUx, _resBig, 0, _resBig[2]); 
+  FLUID_3D::setZeroY(bigUy, _resBig, 0, _resBig[2]); 
+  FLUID_3D::setZeroZ(bigUz, _resBig, 0, _resBig[2]);
 
   // do the MacCormack advection, with substepping if necessary
   for(int substep = 0; substep < totalSubsteps; substep++)
@@ -704,7 +716,7 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel,
   } // substep
   
   // wipe the density borders
-  FLUID_3D::setZeroBorder(_densityBig, _resBig);
+  FLUID_3D::setZeroBorder(_densityBig, _resBig, 0, _resBig[2]);
     
   // reset texture coordinates now in preparation for next timestep
   // Shouldn't do this before generating the noise because then the 
@@ -724,7 +736,9 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel,
   
 
   _totalStepsBig++;
-}
+}*/
+
+//struct
 
 //////////////////////////////////////////////////////////////////////
 // perform the full turbulence algorithm, including OpenMP 
@@ -747,6 +761,7 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa
 
        memset(_tcTemp, 0, sizeof(float)*_totalCellsSm);
 
+
        // prepare textures
        advectTextureCoordinates(dtOrg, xvel,yvel,zvel, tempBig1, tempBig2);
 
@@ -763,25 +778,38 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa
                if (obstacles[x]) highFreqEnergy[x] = 0.f;
 
        Vec3Int ressm(_xResSm, _yResSm, _zResSm);
-       FLUID_3D::setNeumannX(highFreqEnergy, ressm);
-       FLUID_3D::setNeumannY(highFreqEnergy, ressm);
-       FLUID_3D::setNeumannZ(highFreqEnergy, ressm);
+       FLUID_3D::setNeumannX(highFreqEnergy, ressm, 0 , ressm[2]);
+       FLUID_3D::setNeumannY(highFreqEnergy, ressm, 0 , ressm[2]);
+       FLUID_3D::setNeumannZ(highFreqEnergy, ressm, 0 , ressm[2]);
+
+
+   int threadval = 1;
+#if PARALLEL==1
+  threadval = omp_get_max_threads();
+#endif
+
 
   // parallel region setup
-  float maxVelMagThreads[8] = { -1., -1., -1., -1., -1., -1., -1., -1. };
-#if PARALLEL==1 && !_WIN32
+  // Uses omp_get_max_trheads to get number of required cells.
+  float* maxVelMagThreads = new float[threadval];
+
+  for (int i=0; i<threadval; i++) maxVelMagThreads[i] = -1.0f;
+
+#if PARALLEL==1
+
 #pragma omp parallel
 #endif
   { float maxVelMag1 = 0.;
-#if PARALLEL==1 && !_WIN32
+#if PARALLEL==1
     const int id  = omp_get_thread_num(); /*, num = omp_get_num_threads(); */
 #endif
 
   // vector noise main loop
-#if PARALLEL==1 && !_WIN32
-#pragma omp for  schedule(static)
+#if PARALLEL==1
+#pragma omp for schedule(static,1)
 #endif
-  for (int zSmall = 0; zSmall < _zResSm; zSmall++) 
+  for (int zSmall = 0; zSmall < _zResSm; zSmall++)
+  {
   for (int ySmall = 0; ySmall < _yResSm; ySmall++) 
   for (int xSmall = 0; xSmall < _xResSm; xSmall++)
   {
@@ -796,14 +824,14 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa
 
     // get LU factorization of texture jacobian and apply 
     // it to unit vectors
-    JAMA::LU<float> LU = computeLU3x3(jacobian);
+    sLU LU = computeLU(jacobian);
     float xUnwarped[] = {1.0f, 0.0f, 0.0f};
     float yUnwarped[] = {0.0f, 1.0f, 0.0f};
     float zUnwarped[] = {0.0f, 0.0f, 1.0f};
     float xWarped[] = {1.0f, 0.0f, 0.0f};
     float yWarped[] = {0.0f, 1.0f, 0.0f};
     float zWarped[] = {0.0f, 0.0f, 1.0f};
-    bool nonSingular = LU.isNonsingular();
+    bool nonSingular = isNonsingular(LU);
 #if 0
        // UNUSED
     float eigMax = 10.0f;
@@ -908,24 +936,27 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa
           obstacles, posSm[0], posSm[1], posSm[2], _xResSm, _yResSm, _zResSm); 
       if (obsCheck > 0.95) 
         bigUx[index] = bigUy[index] = bigUz[index] = 0.;
-    } // xyz
+    } // xyz*/
 
-#if PARALLEL==1 && !_WIN32
+#if PARALLEL==1
     maxVelMagThreads[id] = maxVelMag1;
 #else
     maxVelMagThreads[0] = maxVelMag1;
 #endif
+  }
   }
   } // omp
   
   // compute maximum over threads
   float maxVelMag = maxVelMagThreads[0];
-#if PARALLEL==1 && !_WIN32
-  for (int i = 1; i < 8; i++) 
+#if PARALLEL==1
+  for (int i = 1; i < threadval; i++) 
     if (maxVelMag < maxVelMagThreads[i]) 
       maxVelMag = maxVelMagThreads[i];
+  delete [] maxVelMagThreads;
 #endif
 
+
   // prepare density for an advection
   SWAP_POINTERS(_densityBig, _densityBigOld);
 
@@ -941,15 +972,56 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa
   const float dtSubdiv = dt / (float)totalSubsteps;
 
   // set boundaries of big velocity grid
-  FLUID_3D::setZeroX(bigUx, _resBig); 
-  FLUID_3D::setZeroY(bigUy, _resBig); 
-  FLUID_3D::setZeroZ(bigUz, _resBig);
+  FLUID_3D::setZeroX(bigUx, _resBig, 0 , _resBig[2]); 
+  FLUID_3D::setZeroY(bigUy, _resBig, 0 , _resBig[2]); 
+  FLUID_3D::setZeroZ(bigUz, _resBig, 0 , _resBig[2]);
+
+#if PARALLEL==1
+  int stepParts = threadval*2; // Dividing parallelized sections into numOfThreads * 2 sections
+  float partSize = (float)_zResBig/stepParts;  // Size of one part;
+
+  if (partSize < 4) {stepParts = threadval;                                    // If the slice gets too low (might actually slow things down, change it to larger
+                                       partSize = (float)_zResBig/stepParts;}
+  if (partSize < 4) {stepParts = (int)(ceil((float)_zResBig/4.0f));    // If it's still too low (only possible on future systems with +24 cores), change it to 4
+                                       partSize = (float)_zResBig/stepParts;}
+#else
+  int zBegin=0;
+  int zEnd=_resBig[2];
+#endif
 
   // do the MacCormack advection, with substepping if necessary
   for(int substep = 0; substep < totalSubsteps; substep++)
   {
-    FLUID_3D::advectFieldMacCormack(dtSubdiv, bigUx, bigUy, bigUz, 
-        _densityBigOld, _densityBig, tempBig1, tempBig2, _resBig, NULL);
+
+#if PARALLEL==1
+       #pragma omp parallel
+       {
+
+       #pragma omp for schedule(static,1)
+       for (int i=0; i<stepParts; i++)
+       {
+               int zBegin = (int)((float)i*partSize + 0.5f);
+               int zEnd = (int)((float)(i+1)*partSize + 0.5f);
+#endif
+               FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, 
+                   _densityBigOld, tempBig1, _resBig, zBegin, zEnd);
+#if PARALLEL==1
+       }
+
+       #pragma omp barrier
+
+       #pragma omp for schedule(static,1)
+       for (int i=0; i<stepParts; i++)
+       {
+               int zBegin = (int)((float)i*partSize + 0.5f);
+               int zEnd = (int)((float)(i+1)*partSize + 0.5f);
+#endif
+               FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, 
+                   _densityBigOld, _densityBig, tempBig1, tempBig2, _resBig, NULL, zBegin, zEnd);
+#if PARALLEL==1
+       }
+       }
+#endif
 
     if (substep < totalSubsteps - 1) 
       SWAP_POINTERS(_densityBig, _densityBigOld);
@@ -964,7 +1036,7 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa
   free(highFreqEnergy);
   
   // wipe the density borders
-  FLUID_3D::setZeroBorder(_densityBig, _resBig);
+  FLUID_3D::setZeroBorder(_densityBig, _resBig, 0 , _resBig[2]);
     
   // reset texture coordinates now in preparation for next timestep
   // Shouldn't do this before generating the noise because then the 
index 2d1d590fcc08cdbabe85ee7a3f411561b7e168fc..c6a3985c8d0b8ed97f8fd20034ba1548fb090b71 100644 (file)
@@ -30,6 +30,7 @@
 
 #include <stdio.h>
 #include <stdlib.h>
+#include <math.h>
 
 // y in smoke is z in blender
 extern "C" FLUID_3D *smoke_init(int *res, float *p0, float dt)
@@ -110,8 +111,8 @@ extern "C" void smoke_dissolve(FLUID_3D *fluid, int speed, int log)
 
                        heat[i] *= (1.0 - dydx);
 
-                       if(heat[i] < 0.0f)
-                               heat[i] = 0.0f;
+                       /*if(heat[i] < 0.0f)
+                               heat[i] = 0.0f;*/
                }
        }
        else // linear falloff
@@ -127,10 +128,9 @@ extern "C" void smoke_dissolve(FLUID_3D *fluid, int speed, int log)
                        if(density[i] < 0.0f)
                                density[i] = 0.0f;
 
-                       heat[i] -= dydx;
-
-                       if(heat[i] < 0.0f)
-                               heat[i] = 0.0f;
+                       if(abs(heat[i]) < dydx) heat[i] = 0.0f;
+                       else if (heat[i]>0.0f) heat[i] -= dydx;
+                       else if (heat[i]<0.0f) heat[i] += dydx;
                                
                }
        }
index 646b1a853dc84e00a9a8df8f49b078b9ab3e3595..8b16452ed34500788b3facb3e6acc6856a2231d1 100644 (file)
@@ -39,6 +39,9 @@
 #include <stdlib.h> 
 #include <ctype.h>
 #include <string.h>
+#if defined(__sun__) || defined( __sun ) || defined (__sparc) || defined (__sparc__) || defined (_AIX)
+#include <strings.h>
+#endif
 #include "STR_String.h"
 
 /*-------------------------------------------------------------------------------------------------
index 68feafe6bef082cff90fe1b778c49693025fdf20..2a68201c8f559c1e6a102521388d403b4845523a 100644 (file)
                                        RelativePath="..\..\..\source\blender\python\intern\bpy_operator_wrap.c"\r
                                        >\r
                                </File>\r
+                               <File\r
+                                       RelativePath="..\..\..\source\blender\python\intern\bpy_props.c"\r
+                                       >\r
+                               </File>\r
                                <File\r
                                        RelativePath="..\..\..\source\blender\python\intern\bpy_rna.c"\r
                                        >\r
                                        RelativePath="..\..\..\source\blender\python\intern\bpy_panel_wrap.h"\r
                                        >\r
                                </File>\r
+                               <File\r
+                                       RelativePath="..\..\..\source\blender\python\intern\bpy_props.h"\r
+                                       >\r
+                               </File>\r
                                <File\r
                                        RelativePath="..\..\..\source\blender\python\intern\bpy_rna.h"\r
                                        >\r
index 7ac0f473b0da38e95da1c8e3bb5f78639555215a..3d6c954a6a9972405be17d5143a9729eaacc8fb9 100644 (file)
                                RelativePath="..\..\..\source\blender\blenlib\intern\math_base.c"\r
                                >\r
                        </File>\r
+                       <File\r
+                               RelativePath="..\..\..\source\blender\blenlib\intern\math_base_inline.c"\r
+                               >\r
+                       </File>\r
                        <File\r
                                RelativePath="..\..\..\source\blender\blenlib\intern\math_color.c"\r
                                >\r
                                RelativePath="..\..\..\source\blender\blenlib\BLI_math_geom.h"\r
                                >\r
                        </File>\r
+                       <File\r
+                               RelativePath="..\..\..\source\blender\blenlib\BLI_math_inline.h"\r
+                               >\r
+                       </File>\r
                        <File\r
                                RelativePath="..\..\..\source\blender\blenlib\BLI_math_matrix.h"\r
                                >\r
index d8e8f5a9eca88196ab7dbb877db2f7130b174082..7d45f96a3b8c61f211924f167339f824fdda183e 100644 (file)
@@ -44,7 +44,7 @@
                                Name="VCCLCompilerTool"\r
                                Optimization="2"\r
                                InlineFunctionExpansion="1"\r
-                               AdditionalIncludeDirectories="..\..\..\..\lib\windows\QTDevWin\CIncludes;..\..\..\..\lib\windows\sdl\include;..\..\..\..\lib\windows\python\include\python2.6;..\..\..\..\build\msvc_9\intern\bsp\include;..\..\..\..\build\msvc_9\intern\ghost\include;..\..\..\..\build\msvc_9\intern\elbeem\include;..\..\..\..\build\msvc_9\intern\opennl\include;..\..\..\..\build\msvc_9\intern\bmfont\include;..\..\..\..\build\msvc_9\intern\blenkey\include;..\..\..\..\build\msvc_9\intern\decimation\include;..\..\..\..\build\msvc_9\intern\memutil\include;..\..\..\..\build\msvc_9\intern\guardedalloc\include;..\..\..\..\build\msvc_9\intern\smoke\include;..\..\..\..\build\msvc_9\intern\soundsystem\include;..\..\..\source\blender;..\..\..\source\blender\img;..\..\..\source\blender\verify;..\..\..\source\blender\ftfont;..\..\..\source\blender\misc;..\..\..\source\blender\imbuf;..\..\..\source\blender\blenlib;..\..\..\source\blender\python;..\..\..\source\blender\editors\include;..\..\..\source\blender\renderui;..\..\..\source\blender\blenloader;..\..\..\source\blender\quicktime;..\..\..\source\blender\blenkernel;..\..\..\source\blender\blenfont;..\..\..\source\blender\makesdna;..\..\..\source\blender\makesrna;..\..\..\source\blender\nodes;..\..\..\source\blender\windowmanager;..\..\..\source\blender\blenpluginapi;..\..\..\source\blender\renderconverter;..\..\..\source\blender\readstreamglue;..\..\..\source\blender\render\extern\include;..\..\..\source\blender\radiosity\extern\include;..\..\..\source\kernel\gen_system;..\..\..\source\gameengine\network;..\..\..\source\gameengine\soundsystem\snd_openal;..\..\..\..\build\msvc_9\extern\verse\include;..\..\..\..\lib\windows\pthreads\include;..\..\..\..\lib\windows\ffmpeg\include;..\..\..\..\build\msvc_9\extern\glew\include;..\..\..\source\blender\gpu;..\..\..\intern\audaspace\intern;..\..\..\source\blender\ikplugin"\r
+                               AdditionalIncludeDirectories="..\..\..\..\lib\windows\QTDevWin\CIncludes;..\..\..\..\lib\windows\sdl\include;..\..\..\..\lib\windows\python\include\python2.6;..\..\..\..\build\msvc_9\intern\bsp\include;..\..\..\..\build\msvc_9\intern\ghost\include;..\..\..\..\build\msvc_9\intern\elbeem\include;..\..\..\..\build\msvc_9\intern\opennl\include;..\..\..\..\build\msvc_9\intern\bmfont\include;..\..\..\..\build\msvc_9\intern\blenkey\include;..\..\..\..\build\msvc_9\intern\decimation\include;..\..\..\..\build\msvc_9\intern\memutil\include;..\..\..\..\build\msvc_9\intern\guardedalloc\include;..\..\..\..\build\msvc_9\intern\smoke\include;..\..\..\..\build\msvc_9\intern\soundsystem\include;..\..\..\source\blender;..\..\..\source\blender\img;..\..\..\source\blender\verify;..\..\..\source\blender\ftfont;..\..\..\source\blender\misc;..\..\..\source\blender\imbuf;..\..\..\source\blender\blenlib;..\..\..\source\blender\python;..\..\..\source\blender\editors\include;..\..\..\source\blender\editors\interface;..\..\..\source\blender\renderui;..\..\..\source\blender\blenloader;..\..\..\source\blender\quicktime;..\..\..\source\blender\blenkernel;..\..\..\source\blender\blenfont;..\..\..\source\blender\makesdna;..\..\..\source\blender\makesrna;..\..\..\source\blender\nodes;..\..\..\source\blender\windowmanager;..\..\..\source\blender\blenpluginapi;..\..\..\source\blender\renderconverter;..\..\..\source\blender\readstreamglue;..\..\..\source\blender\render\extern\include;..\..\..\source\blender\radiosity\extern\include;..\..\..\source\kernel\gen_system;..\..\..\source\gameengine\network;..\..\..\source\gameengine\soundsystem\snd_openal;..\..\..\..\build\msvc_9\extern\verse\include;..\..\..\..\lib\windows\pthreads\include;..\..\..\..\lib\windows\ffmpeg\include;..\..\..\..\build\msvc_9\extern\glew\include;..\..\..\source\blender\gpu;..\..\..\intern\audaspace\intern;..\..\..\source\blender\ikplugin"\r
                                PreprocessorDefinitions="NDEBUG;WIN32;_LIB;_CONSOLE;GAMEBLENDER=1;WITH_QUICKTIME;INTERNATIONAL;WITH_FREETYPE2;WITH_INTERNATIONAL;WITH_OPENEXR;WITH_DDS;WITH_BULLET=1;WITH_FFMPEG"\r
                                StringPooling="true"\r
                                RuntimeLibrary="0"\r
                        <Tool\r
                                Name="VCCLCompilerTool"\r
                                Optimization="0"\r
-                               AdditionalIncludeDirectories="..\..\..\..\lib\windows\QTDevWin\CIncludes;..\..\..\..\lib\windows\sdl\include;..\..\..\..\lib\windows\python\include\python2.6;..\..\..\..\build\msvc_9\intern\bsp\include;..\..\..\..\build\msvc_9\intern\ghost\include;..\..\..\..\build\msvc_9\intern\elbeem\include;..\..\..\..\build\msvc_9\intern\opennl\include;..\..\..\..\build\msvc_9\intern\bmfont\include;..\..\..\..\build\msvc_9\intern\blenkey\include;..\..\..\..\build\msvc_9\intern\decimation\include;..\..\..\..\build\msvc_9\intern\memutil\include;..\..\..\..\build\msvc_9\intern\guardedalloc\include;..\..\..\..\build\msvc_9\intern\smoke\include;..\..\..\..\build\msvc_9\intern\soundsystem\include;..\..\..\source\blender;..\..\..\source\blender\img;..\..\..\source\blender\verify;..\..\..\source\blender\ftfont;..\..\..\source\blender\misc;..\..\..\source\blender\imbuf;..\..\..\source\blender\blenlib;..\..\..\source\blender\python;..\..\..\source\blender\editors\include;..\..\..\source\blender\renderui;..\..\..\source\blender\blenloader;..\..\..\source\blender\quicktime;..\..\..\source\blender\blenkernel;..\..\..\source\blender\blenfont;..\..\..\source\blender\makesdna;..\..\..\source\blender\makesrna;..\..\..\source\blender\nodes;..\..\..\source\blender\windowmanager;..\..\..\source\blender\blenpluginapi;..\..\..\source\blender\renderconverter;..\..\..\source\blender\readstreamglue;..\..\..\source\blender\render\extern\include;..\..\..\source\blender\radiosity\extern\include;..\..\..\source\kernel\gen_system;..\..\..\source\gameengine\network;..\..\..\source\gameengine\soundsystem\snd_openal;..\..\..\..\build\msvc_9\extern\verse\include;..\..\..\..\lib\windows\pthreads\include;..\..\..\..\lib\windows\ffmpeg\include;..\..\..\..\build\msvc_9\extern\glew\include;..\..\..\source\blender\gpu;..\..\..\intern\audaspace\intern;..\..\..\source\blender\ikplugin"\r
+                               AdditionalIncludeDirectories="..\..\..\..\lib\windows\QTDevWin\CIncludes;..\..\..\..\lib\windows\sdl\include;..\..\..\..\lib\windows\python\include\python2.6;..\..\..\..\build\msvc_9\intern\bsp\include;..\..\..\..\build\msvc_9\intern\ghost\include;..\..\..\..\build\msvc_9\intern\elbeem\include;..\..\..\..\build\msvc_9\intern\opennl\include;..\..\..\..\build\msvc_9\intern\bmfont\include;..\..\..\..\build\msvc_9\intern\blenkey\include;..\..\..\..\build\msvc_9\intern\decimation\include;..\..\..\..\build\msvc_9\intern\memutil\include;..\..\..\..\build\msvc_9\intern\guardedalloc\include;..\..\..\..\build\msvc_9\intern\smoke\include;..\..\..\..\build\msvc_9\intern\soundsystem\include;..\..\..\source\blender;..\..\..\source\blender\img;..\..\..\source\blender\verify;..\..\..\source\blender\ftfont;..\..\..\source\blender\misc;..\..\..\source\blender\imbuf;..\..\..\source\blender\blenlib;..\..\..\source\blender\python;..\..\..\source\blender\editors\include;..\..\..\source\blender\editors\interface;..\..\..\source\blender\renderui;..\..\..\source\blender\blenloader;..\..\..\source\blender\quicktime;..\..\..\source\blender\blenkernel;..\..\..\source\blender\blenfont;..\..\..\source\blender\makesdna;..\..\..\source\blender\makesrna;..\..\..\source\blender\nodes;..\..\..\source\blender\windowmanager;..\..\..\source\blender\blenpluginapi;..\..\..\source\blender\renderconverter;..\..\..\source\blender\readstreamglue;..\..\..\source\blender\render\extern\include;..\..\..\source\blender\radiosity\extern\include;..\..\..\source\kernel\gen_system;..\..\..\source\gameengine\network;..\..\..\source\gameengine\soundsystem\snd_openal;..\..\..\..\build\msvc_9\extern\verse\include;..\..\..\..\lib\windows\pthreads\include;..\..\..\..\lib\windows\ffmpeg\include;..\..\..\..\build\msvc_9\extern\glew\include;..\..\..\source\blender\gpu;..\..\..\intern\audaspace\intern;..\..\..\source\blender\ikplugin"\r
                                PreprocessorDefinitions="_DEBUG;WIN32;_LIB;_CONSOLE;GAMEBLENDER=1;WITH_QUICKTIME;INTERNATIONAL;WITH_BF_INTERNATIONAL;WITH_FREETYPE2;WITH_OPENEXR;WITH_DDS;WITH_BULLET = 1;WITH_FFMPEG"\r
                                BasicRuntimeChecks="3"\r
                                RuntimeLibrary="1"\r
                <Filter\r
                        Name="space_view3d"\r
                        >\r
+                       <File\r
+                               RelativePath="..\..\..\source\blender\editors\space_view3d\drawanimviz.c"\r
+                               >\r
+                       </File>\r
                        <File\r
                                RelativePath="..\..\..\source\blender\editors\space_view3d\drawarmature.c"\r
                                >\r
index e9cb51b5cff62aceba4ee817eab6d3489ef5ef1a..a86e91670c49ebef745801608bfdde6d6ba0a0d6 100644 (file)
                                        RelativePath="..\..\..\source\blender\nodes\intern\CMP_nodes\CMP_chromaMatte.c"\r
                                        >\r
                                </File>\r
+                               <File\r
+                                       RelativePath="..\..\..\source\blender\nodes\intern\CMP_nodes\CMP_colorbalance.c"\r
+                                       >\r
+                               </File>\r
                                <File\r
                                        RelativePath="..\..\..\source\blender\nodes\intern\CMP_nodes\CMP_colorMatte.c"\r
                                        >\r
                                        RelativePath="..\..\..\source\blender\nodes\intern\CMP_nodes\CMP_glare.c"\r
                                        >\r
                                </File>\r
+                               <File\r
+                                       RelativePath="..\..\..\source\blender\nodes\intern\CMP_nodes\CMP_huecorrect.c"\r
+                                       >\r
+                               </File>\r
                                <File\r
                                        RelativePath="..\..\..\source\blender\nodes\intern\CMP_nodes\CMP_hueSatVal.c"\r
                                        >\r
index af046f25563fe0b17f7352fe527ad0ab2058f049..2aee30b407fcc0b7cfb64d8abf67ec4090cccf7d 100644 (file)
Binary files a/release/datafiles/blenderbuttons and b/release/datafiles/blenderbuttons differ
index 5275574d903930a26146f3dce85f78c02f0e6b7e..7abb0c19cd7debf236b55e88e216a85977e91957 100644 (file)
@@ -146,7 +146,7 @@ def write_pov(filename, scene=None, info_callback=None):
         file.write('\tup <0, 1, 0>\n')
         file.write('\tangle  %f \n' % (360.0 * atan(16.0 / camera.data.lens) / pi))
 
-        file.write('\trotate  <%.6f, %.6f, %.6f>\n' % tuple([degrees(e) for e in matrix.rotationPart().toEuler()]))
+        file.write('\trotate  <%.6f, %.6f, %.6f>\n' % tuple([degrees(e) for e in matrix.rotation_part().to_euler()]))
         file.write('\ttranslate <%.6f, %.6f, %.6f>\n' % (matrix[3][0], matrix[3][1], matrix[3][2]))
         file.write('}\n')
 
index 979c53c581d744bf1fb0398ee7f786ebcbc4665c..f898a4b5f9ea344121d7f873763a2189efa31881 100644 (file)
@@ -839,7 +839,7 @@ def make_track_chunk(ID, obj):
             track_chunk.add_variable("position", _3ds_point_3d(obj.getLocation()))
         elif ID==ROT_TRACK_TAG:
             # rotation (quaternion, angle first, followed by axis):
-            q = obj.getEuler().toQuat()
+            q = obj.getEuler().to_quat()
             track_chunk.add_variable("rotation", _3ds_point_4d((q.angle, q.axis[0], q.axis[1], q.axis[2])))
         elif ID==SCL_TRACK_TAG:
             # scale vector:
index ffd7c3f3efa6f6141eac8cc366d53931bd5b8604..dcf109e24e1e8e8a7ce0920e15c717df916a792c 100644 (file)
@@ -146,7 +146,7 @@ def eulerRadToDeg(eul):
 mtx4_identity = Mathutils.Matrix()
 
 # testing
-mtx_x90                = Mathutils.RotationMatrix( math.pi/2, 3, 'x') # used
+mtx_x90                = Mathutils.RotationMatrix( math.pi/2, 3, 'X') # used
 #mtx_x90n      = RotationMatrix(-90, 3, 'x')
 #mtx_y90       = RotationMatrix( 90, 3, 'y')
 #mtx_y90n      = RotationMatrix(-90, 3, 'y')
@@ -154,11 +154,11 @@ mtx_x90           = Mathutils.RotationMatrix( math.pi/2, 3, 'x') # used
 #mtx_z90n      = RotationMatrix(-90, 3, 'z')
 
 #mtx4_x90      = RotationMatrix( 90, 4, 'x')
-mtx4_x90n      = Mathutils.RotationMatrix(-math.pi/2, 4, 'x') # used
+mtx4_x90n      = Mathutils.RotationMatrix(-math.pi/2, 4, 'X') # used
 #mtx4_y90      = RotationMatrix( 90, 4, 'y')
-mtx4_y90n      = Mathutils.RotationMatrix(-math.pi/2, 4, 'y') # used
-mtx4_z90       = Mathutils.RotationMatrix( math.pi/2, 4, 'z') # used
-mtx4_z90n      = Mathutils.RotationMatrix(-math.pi/2, 4, 'z') # used
+mtx4_y90n      = Mathutils.RotationMatrix(-math.pi/2, 4, 'Y') # used
+mtx4_z90       = Mathutils.RotationMatrix( math.pi/2, 4, 'Z') # used
+mtx4_z90n      = Mathutils.RotationMatrix(-math.pi/2, 4, 'Z') # used
 
 # def strip_path(p):
 #      return p.split('\\')[-1].split('/')[-1]
@@ -590,9 +590,9 @@ def write(filename, batch_objects = None, \
         def getAnimParRelMatrixRot(self, frame):
             type = self.blenObject.type
             if self.fbxParent:
-                matrix_rot = (((self.__anim_poselist[frame] * GLOBAL_MATRIX) * (self.fbxParent.__anim_poselist[frame] * GLOBAL_MATRIX).invert())).rotationPart()
+                matrix_rot = (((self.__anim_poselist[frame] * GLOBAL_MATRIX) * (self.fbxParent.__anim_poselist[frame] * GLOBAL_MATRIX).invert())).rotation_part()
             else:
-                matrix_rot = (self.__anim_poselist[frame] * GLOBAL_MATRIX).rotationPart()
+                matrix_rot = (self.__anim_poselist[frame] * GLOBAL_MATRIX).rotation_part()
 
             # Lamps need to be rotated
             if type =='LAMP':
@@ -600,7 +600,7 @@ def write(filename, batch_objects = None, \
             elif type =='CAMERA':
 #                      elif ob and type =='Camera':
                 y = Mathutils.Vector(0,1,0) * matrix_rot
-                matrix_rot = matrix_rot * Mathutils.RotationMatrix(math.pi/2, 3, 'r', y)
+                matrix_rot = matrix_rot * Mathutils.RotationMatrix(math.pi/2, 3, y)
 
             return matrix_rot
 
@@ -676,11 +676,11 @@ def write(filename, batch_objects = None, \
 #                              par_matrix = mtx4_z90 * parent.matrix['ARMATURESPACE'] # dont apply armature matrix anymore
                 matrix = matrix * par_matrix.copy().invert()
 
-            matrix_rot =       matrix.rotationPart()
+            matrix_rot =       matrix.rotation_part()
 
-            loc =                      tuple(matrix.translationPart())
-            scale =                    tuple(matrix.scalePart())
-            rot =                      tuple(matrix_rot.toEuler())
+            loc =                      tuple(matrix.translation_part())
+            scale =                    tuple(matrix.scale_part())
+            rot =                      tuple(matrix_rot.to_euler())
 
         else:
             # This is bad because we need the parent relative matrix from the fbx parent (if we have one), dont use anymore
@@ -692,20 +692,20 @@ def write(filename, batch_objects = None, \
             #  matrix = matrix_scale * matrix
 
             if matrix:
-                loc = tuple(matrix.translationPart())
-                scale = tuple(matrix.scalePart())
+                loc = tuple(matrix.translation_part())
+                scale = tuple(matrix.scale_part())
 
-                matrix_rot = matrix.rotationPart()
+                matrix_rot = matrix.rotation_part()
                 # Lamps need to be rotated
                 if ob and ob.type =='Lamp':
                     matrix_rot = mtx_x90 * matrix_rot
-                    rot = tuple(matrix_rot.toEuler())
+                    rot = tuple(matrix_rot.to_euler())
                 elif ob and ob.type =='Camera':
                     y = Mathutils.Vector(0,1,0) * matrix_rot
-                    matrix_rot = matrix_rot * Mathutils.RotationMatrix(math.pi/2, 3, 'r', y)
-                    rot = tuple(matrix_rot.toEuler())
+                    matrix_rot = matrix_rot * Mathutils.RotationMatrix(math.pi/2, 3, y)
+                    rot = tuple(matrix_rot.to_euler())
                 else:
-                    rot = tuple(matrix_rot.toEuler())
+                    rot = tuple(matrix_rot.to_euler())
             else:
                 if not loc:
                     loc = 0,0,0
@@ -1131,7 +1131,7 @@ def write(filename, batch_objects = None, \
         else:
             do_light = 1
 
-        scale = abs(GLOBAL_MATRIX.scalePart()[0]) # scale is always uniform in this case
+        scale = abs(GLOBAL_MATRIX.scale_part()[0]) # scale is always uniform in this case
 
         file.write('\n\t\t\tProperty: "LightType", "enum", "",%i' % light_type)
         file.write('\n\t\t\tProperty: "CastLightOnObject", "bool", "",1')
@@ -1753,7 +1753,7 @@ def write(filename, batch_objects = None, \
                 for uf in uvlayer.data:
 #               for f in me.faces:
                     # workaround, since uf.uv iteration is wrong atm
-                    for uv in [uf.uv1, uf.uv2, uf.uv3, uf.uv4][:len(uf.uv)]:
+                    for uv in uf.uv:
 #                   for uv in f.uv:
                         if i==-1:
                             file.write('%.6f,%.6f' % tuple(uv))
@@ -2866,18 +2866,18 @@ Takes:  {''')
                         # ----------------
                         for TX_LAYER, TX_CHAN in enumerate('TRS'): # transform, rotate, scale
 
-                            if         TX_CHAN=='T':   context_bone_anim_vecs = [mtx[0].translationPart()      for mtx in context_bone_anim_mats]
-                            elif       TX_CHAN=='S':   context_bone_anim_vecs = [mtx[0].scalePart()            for mtx in context_bone_anim_mats]
+                            if         TX_CHAN=='T':   context_bone_anim_vecs = [mtx[0].translation_part()     for mtx in context_bone_anim_mats]
+                            elif       TX_CHAN=='S':   context_bone_anim_vecs = [mtx[0].scale_part()           for mtx in context_bone_anim_mats]
                             elif       TX_CHAN=='R':
                                 # Was....
-                                # elif         TX_CHAN=='R':   context_bone_anim_vecs = [mtx[1].toEuler()                      for mtx in context_bone_anim_mats]
+                                # elif         TX_CHAN=='R':   context_bone_anim_vecs = [mtx[1].to_euler()                     for mtx in context_bone_anim_mats]
                                 #
                                 # ...but we need to use the previous euler for compatible conversion.
                                 context_bone_anim_vecs = []
                                 prev_eul = None
                                 for mtx in context_bone_anim_mats:
-                                    if prev_eul:       prev_eul = mtx[1].toEuler(prev_eul)
-                                    else:                      prev_eul = mtx[1].toEuler()
+                                    if prev_eul:       prev_eul = mtx[1].to_euler(prev_eul)
+                                    else:                      prev_eul = mtx[1].to_euler()
                                     context_bone_anim_vecs.append(eulerRadToDeg(prev_eul))
 #                                                                      context_bone_anim_vecs.append(prev_eul)
 
index 7c4c83c53c7132f846de6ff92646d0bea8b4f7e0..4f3212fa72531d01800e0ede783b0bbf486d0fb6 100644 (file)
@@ -564,7 +564,7 @@ def write(filename, objects, scene,
                     tface = uv_layer.data[f_index]
 
                     # workaround, since tface.uv iteration is wrong atm
-                    uvs = [tface.uv1, tface.uv2, tface.uv3, tface.uv4][:len(tface.uv)]
+                    uvs = tface.uv
                     # uvs = [tface.uv1, tface.uv2, tface.uv3]
 
                     # # add another UV if it's a quad
index c492c4a6b1558ee1f886659f40c98362b6d6fc6e..0875abf1de1cec37baaf6ecd1fb349a6b0c88247 100644 (file)
@@ -81,7 +81,7 @@ from export_3ds import create_derived_objects, free_derived_objects
 
 #
 DEG2RAD=0.017453292519943295
-MATWORLD= Mathutils.RotationMatrix(-90, 4, 'x')
+MATWORLD= Mathutils.RotationMatrix(-90, 4, 'X')
 
 ####################################
 # Global Variables
@@ -239,8 +239,8 @@ class x3d_class:
         # get the camera location, subtract 90 degress from X to orient like X3D does
         # mat = ob.matrixWorld - mat is now passed!
 
-        loc = self.rotatePointForVRML(mat.translationPart())
-        rot = mat.toEuler()
+        loc = self.rotatePointForVRML(mat.translation_part())
+        rot = mat.to_euler()
         rot = (((rot[0]-90)), rot[1], rot[2])
         # rot = (((rot[0]-90)*DEG2RAD), rot[1]*DEG2RAD, rot[2]*DEG2RAD)
         nRot = self.rotatePointForVRML( rot )
@@ -300,8 +300,8 @@ class x3d_class:
         # note -dz seems to equal om[3][1]
         # note  dy seems to equal om[3][2]
 
-        #location=(ob.matrixWorld*MATWORLD).translationPart() # now passed
-        location=(mtx*MATWORLD).translationPart()
+        #location=(ob.matrixWorld*MATWORLD).translation_part() # now passed
+        location=(mtx*MATWORLD).translation_part()
 
         radius = lamp.distance*math.cos(beamWidth)
         # radius = lamp.dist*math.cos(beamWidth)
@@ -346,8 +346,8 @@ class x3d_class:
             ambi = 0
             ambientIntensity = 0
 
-        # location=(ob.matrixWorld*MATWORLD).translationPart() # now passed
-        location= (mtx*MATWORLD).translationPart()
+        # location=(ob.matrixWorld*MATWORLD).translation_part() # now passed
+        location= (mtx*MATWORLD).translation_part()
 
         self.file.write("<PointLight DEF=\"%s\" " % safeName)
         self.file.write("ambientIntensity=\"%s\" " % (round(ambientIntensity,self.cp)))
@@ -364,8 +364,8 @@ class x3d_class:
             return
         else:
             dx,dy,dz = self.computeDirection(mtx)
-            # location=(ob.matrixWorld*MATWORLD).translationPart()
-            location=(mtx*MATWORLD).translationPart()
+            # location=(ob.matrixWorld*MATWORLD).translation_part()
+            location=(mtx*MATWORLD).translation_part()
             self.writeIndented("<%s\n" % obname,1)
             self.writeIndented("direction=\"%s %s %s\"\n" % (round(dx,3),round(dy,3),round(dz,3)))
             self.writeIndented("location=\"%s %s %s\"\n" % (round(location[0],3), round(location[1],3), round(location[2],3)))
@@ -448,9 +448,9 @@ class x3d_class:
         # mtx = ob.matrixWorld * MATWORLD # mtx is now passed
         mtx = mtx * MATWORLD
 
-        loc= mtx.translationPart()
-        sca= mtx.scalePart()
-        quat = mtx.toQuat()
+        loc= mtx.translation_part()
+        sca= mtx.scale_part()
+        quat = mtx.to_quat()
         rot= quat.axis
 
         self.writeIndented('<Transform DEF="%s" translation="%.6f %.6f %.6f" scale="%.6f %.6f %.6f" rotation="%.6f %.6f %.6f %.6f">\n' % \
@@ -617,7 +617,7 @@ class x3d_class:
         for face in mesh.active_uv_texture.data:
         # for face in mesh.faces:
             # workaround, since tface.uv iteration is wrong atm
-            uvs = [face.uv1, face.uv2, face.uv3, face.uv4][:len(face.uv)]
+            uvs = face.uv
             # uvs = [face.uv1, face.uv2, face.uv3, face.uv4] if face.verts[3] else [face.uv1, face.uv2, face.uv3]
 
             for uv in uvs:
@@ -1048,7 +1048,7 @@ class x3d_class:
     def computeDirection(self, mtx):
         x,y,z=(0,-1.0,0) # point down
 
-        ax,ay,az = (mtx*MATWORLD).toEuler()
+        ax,ay,az = (mtx*MATWORLD).to_euler()
 
         # ax *= DEG2RAD
         # ay *= DEG2RAD
index 2afd2db3cb2e847af3db03315aef4b89691b75f3..f6626e8fb1eef3a0e91f3868cfde4ef2a8150e88 100644 (file)
@@ -89,13 +89,13 @@ MATRIX_IDENTITY_4x4 = Matrix([1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1])
 def eulerRotate(x,y,z, rot_order):
 
     # Clamp all values between 0 and 360, values outside this raise an error.
-    mats=[RotationMatrix(math.radians(x%360),3,'x'), RotationMatrix(math.radians(y%360),3,'y'), RotationMatrix(math.radians(z%360),3,'z')]
+    mats=[RotationMatrix(math.radians(x % 360), 3, 'X'), RotationMatrix(math.radians(y % 360),3,'Y'), RotationMatrix(math.radians(z % 360), 3, 'Z')]
     # print rot_order
     # Standard BVH multiplication order, apply the rotation in the order Z,X,Y
 
     #XXX, order changes???
-    #eul = (mats[rot_order[2]]*(mats[rot_order[1]]* (mats[rot_order[0]]* MATRIX_IDENTITY_3x3))).toEuler()
-    eul = (MATRIX_IDENTITY_3x3*mats[rot_order[0]]*(mats[rot_order[1]]* (mats[rot_order[2]]))).toEuler()
+    #eul = (mats[rot_order[2]]*(mats[rot_order[1]]* (mats[rot_order[0]]* MATRIX_IDENTITY_3x3))).to_euler()
+    eul = (MATRIX_IDENTITY_3x3*mats[rot_order[0]]*(mats[rot_order[1]]* (mats[rot_order[2]]))).to_euler()
 
     eul = math.degrees(eul.x), math.degrees(eul.y), math.degrees(eul.z)
 
@@ -534,8 +534,8 @@ def bvh_node_dict2armature(context, bvh_nodes, IMPORT_START_FRAME= 1, IMPORT_LOO
         bone_name= bvh_node.temp # may not be the same name as the bvh_node, could have been shortened.
         pose_bone= pose_bones[bone_name]
         rest_bone= arm_data.bones[bone_name]
-#XXX           bone_rest_matrix = rest_bone.matrix['ARMATURESPACE'].rotationPart()
-        bone_rest_matrix = rest_bone.matrix.rotationPart()
+#XXX           bone_rest_matrix = rest_bone.matrix['ARMATURESPACE'].rotation_part()
+        bone_rest_matrix = rest_bone.matrix.rotation_part()
 
 
         bone_rest_matrix_inv= Matrix(bone_rest_matrix)
@@ -575,31 +575,31 @@ def bvh_node_dict2armature(context, bvh_nodes, IMPORT_START_FRAME= 1, IMPORT_LOO
 
                 if ROT_STYLE=='QUAT':
                     # Set the rotation, not so simple
-                    bone_rotation_matrix= Euler(math.radians(rx), math.radians(ry), math.radians(rz)).toMatrix()
+                    bone_rotation_matrix= Euler(math.radians(rx), math.radians(ry), math.radians(rz)).to_matrix()
 
                     bone_rotation_matrix.resize4x4()
                     #XXX ORDER CHANGE???
-                    #pose_bone.rotation_quaternion= (bone_rest_matrix * bone_rotation_matrix * bone_rest_matrix_inv).toQuat() # ORIGINAL
-                    # pose_bone.rotation_quaternion= (bone_rest_matrix_inv * bone_rotation_matrix * bone_rest_matrix).toQuat()
-                    # pose_bone.rotation_quaternion= (bone_rotation_matrix * bone_rest_matrix).toQuat() # BAD
-                    # pose_bone.rotation_quaternion= bone_rotation_matrix.toQuat() # NOT GOOD
-                    # pose_bone.rotation_quaternion= bone_rotation_matrix.toQuat() # NOT GOOD
+                    #pose_bone.rotation_quaternion= (bone_rest_matrix * bone_rotation_matrix * bone_rest_matrix_inv).to_quat() # ORIGINAL
+                    # pose_bone.rotation_quaternion= (bone_rest_matrix_inv * bone_rotation_matrix * bone_rest_matrix).to_quat()
+                    # pose_bone.rotation_quaternion= (bone_rotation_matrix * bone_rest_matrix).to_quat() # BAD
+                    # pose_bone.rotation_quaternion= bone_rotation_matrix.to_quat() # NOT GOOD
+                    # pose_bone.rotation_quaternion= bone_rotation_matrix.to_quat() # NOT GOOD
 
-                    #pose_bone.rotation_quaternion= (bone_rotation_matrix * bone_rest_matrix_inv * bone_rest_matrix).toQuat()
-                    #pose_bone.rotation_quaternion= (bone_rest_matrix_inv * bone_rest_matrix * bone_rotation_matrix).toQuat()
-                    #pose_bone.rotation_quaternion= (bone_rest_matrix * bone_rotation_matrix * bone_rest_matrix_inv).toQuat()
+                    #pose_bone.rotation_quaternion= (bone_rotation_matrix * bone_rest_matrix_inv * bone_rest_matrix).to_quat()
+                    #pose_bone.rotation_quaternion= (bone_rest_matrix_inv * bone_rest_matrix * bone_rotation_matrix).to_quat()
+                    #pose_bone.rotation_quaternion= (bone_rest_matrix * bone_rotation_matrix * bone_rest_matrix_inv).to_quat()
 
-                    #pose_bone.rotation_quaternion= ( bone_rest_matrix* bone_rest_matrix_inv * bone_rotation_matrix).toQuat()
-                    #pose_bone.rotation_quaternion= (bone_rotation_matrix * bone_rest_matrix  * bone_rest_matrix_inv).toQuat()
-                    #pose_bone.rotation_quaternion= (bone_rest_matrix_inv * bone_rotation_matrix  * bone_rest_matrix ).toQuat()
+                    #pose_bone.rotation_quaternion= ( bone_rest_matrix* bone_rest_matrix_inv * bone_rotation_matrix).to_quat()
+                    #pose_bone.rotation_quaternion= (bone_rotation_matrix * bone_rest_matrix  * bone_rest_matrix_inv).to_quat()
+                    #pose_bone.rotation_quaternion= (bone_rest_matrix_inv * bone_rotation_matrix  * bone_rest_matrix ).to_quat()
 
-                    pose_bone.rotation_quaternion= (bone_rest_matrix_inv * bone_rotation_matrix * bone_rest_matrix).toQuat()
+                    pose_bone.rotation_quaternion= (bone_rest_matrix_inv * bone_rotation_matrix * bone_rest_matrix).to_quat()
 
                 else:
-                    bone_rotation_matrix= Euler(math.radians(rx), math.radians(ry), math.radians(rz)).toMatrix()
+                    bone_rotation_matrix= Euler(math.radians(rx), math.radians(ry), math.radians(rz)).to_matrix()
                     bone_rotation_matrix.resize4x4()
 
-                    eul= (bone_rest_matrix * bone_rotation_matrix * bone_rest_matrix_inv).toEuler()
+                    eul= (bone_rest_matrix * bone_rotation_matrix * bone_rest_matrix_inv).to_euler()
 
                     #pose_bone.rotation_euler = math.radians(rx), math.radians(ry), math.radians(rz)
                     pose_bone.rotation_euler = eul
@@ -610,8 +610,8 @@ def bvh_node_dict2armature(context, bvh_nodes, IMPORT_START_FRAME= 1, IMPORT_LOO
                 # Set the Location, simple too
 
                 #XXX ORDER CHANGE
-                # pose_bone.location= (TranslationMatrix(Vector(lx, ly, lz) - bvh_node.rest_head_local ) * bone_rest_matrix_inv).translationPart() # WHY * 10? - just how pose works
-                # pose_bone.location= (bone_rest_matrix_inv * TranslationMatrix(Vector(lx, ly, lz) - bvh_node.rest_head_local )).translationPart()
+                # pose_bone.location= (TranslationMatrix(Vector(lx, ly, lz) - bvh_node.rest_head_local ) * bone_rest_matrix_inv).translation_part() # WHY * 10? - just how pose works
+                # pose_bone.location= (bone_rest_matrix_inv * TranslationMatrix(Vector(lx, ly, lz) - bvh_node.rest_head_local )).translation_part()
                 # pose_bone.location= lx, ly, lz
                 pose_bone.location= Vector(lx, ly, lz) - bvh_node.rest_head_local
 
@@ -714,12 +714,12 @@ def bvh_node_dict2armature(context, bvh_nodes, IMPORT_START_FRAME= 1, IMPORT_LOO
 
 
         def pose_rot(anim_data):
-            bone_rotation_matrix= Euler(anim_data[3], anim_data[4], anim_data[5]).toMatrix()
+            bone_rotation_matrix= Euler(anim_data[3], anim_data[4], anim_data[5]).to_matrix()
             bone_rotation_matrix.resize4x4()
-            return tuple((bone_rest_matrix * bone_rotation_matrix * bone_rest_matrix_inv).toQuat()) # qw,qx,qy,qz
+            return tuple((bone_rest_matrix * bone_rotation_matrix * bone_rest_matrix_inv).to_quat()) # qw,qx,qy,qz
 
         def pose_loc(anim_data):
-            return tuple((TranslationMatrix(Vector(anim_data[0], anim_data[1], anim_data[2])) * bone_rest_matrix_inv).translationPart())
+            return tuple((TranslationMatrix(Vector(anim_data[0], anim_data[1], anim_data[2])) * bone_rest_matrix_inv).translation_part())
 
 
         last_frame= len(bvh_node.anim_data)+IMPORT_START_FRAME-1
index 28f082c9bd1fd782b90acf7db3b79934aa9fcb05..a6bb2921a3406fd8cf82dfcf27453c86ae3f3b2c 100644 (file)
@@ -34,8 +34,8 @@ from bpy import ops as _ops_module
 # fake operator module
 ops = _ops_module.ops_fake_module
 
-import sys
-DEBUG = ("-d" in sys.argv)
+import sys as _sys
+DEBUG = ("-d" in _sys.argv)
 
 
 def load_scripts(reload_scripts=False):
@@ -65,6 +65,12 @@ def load_scripts(reload_scripts=False):
             traceback.print_exc()
             return None
 
+    def test_reload(module):
+        try:
+            reload(module)
+        except:
+            traceback.print_exc()
+
     if reload_scripts:
         # reload modules that may not be directly included
         for type_class_name in dir(types):
@@ -76,14 +82,14 @@ def load_scripts(reload_scripts=False):
 
         for module_name in loaded_modules:
             print("Reloading:", module_name)
-            reload(sys.modules[module_name])
+            test_reload(_sys.modules[module_name])
 
     for base_path in utils.script_paths():
         for path_subdir in ("ui", "op", "io"):
             path = os.path.join(base_path, path_subdir)
             if os.path.isdir(path):
-                if path not in sys.path: # reloading would add twice
-                    sys.path.insert(0, path)
+                if path not in _sys.path: # reloading would add twice
+                    _sys.path.insert(0, path)
                 for f in sorted(os.listdir(path)):
                     if f.endswith(".py"):
                         # python module
@@ -96,7 +102,7 @@ def load_scripts(reload_scripts=False):
 
                     if reload_scripts and mod:
                         print("Reloading:", mod)
-                        reload(mod)
+                        test_reload(mod)
 
     if DEBUG:
         print("Time %.4f" % (time.time() - t_main))
@@ -104,10 +110,13 @@ def load_scripts(reload_scripts=False):
 
 def _main():
 
+    # security issue, dont allow the $CWD in the path.
+    _sys.path[:] = filter(None, _sys.path)
+
     # a bit nasty but this prevents help() and input() from locking blender
     # Ideally we could have some way for the console to replace sys.stdin but
     # python would lock blender while waiting for a return value, not easy :|
-    sys.stdin = None
+    _sys.stdin = None
 
     # if "-d" in sys.argv: # Enable this to measure startup speed
     if 0:
index da45ab8eb302876c97952b79a7e5dd760325d02a..8fc92175d0eccf12abe05571cecac357cfb63b74 100644 (file)
 
 # <pep8 compliant>
 
+"""
+This module contains application values that remain unchanged during runtime.
+
+.. data:: version
+
+   The Blender version as a tuple of 3 numbers. eg. (2, 50, 11)
+
+
+.. data:: version_string
+
+   The Blender version formatted as a string.
+
+.. data:: home
+
+   The blender home directory, normally matching $HOME
+
+.. data:: binary_path
+
+   The location of blenders executable, useful for utilities that spawn new instances.
+
+"""
 # constants
 import _bpy
 version = _bpy._VERSION
index 75976ebf9a1834c9b525c813d42593829ecddea6..971a14208a6b67b3ea7dcb68f9ee4534e70b71a9 100644 (file)
@@ -134,7 +134,7 @@ class bpy_ops_submodule_op(object):
 
     def idname(self):
         # submod.foo -> SUBMOD_OT_foo
-        return self.module.upper() + '_OT_' + self.func
+        return self.module + '.' + self.func
 
     def __call__(self, *args, **kw):
 
index 1948a28a7260d13480419e667201a740049bfbea..855751a07438526f389129943bc7c9bcfaf59104 100644 (file)
 
 # <pep8 compliant>
 
-import bpy
-import os
+"""
+This module contains utility functions spesific to blender but
+not assosiated with blenders internal data.
+"""
+
+import bpy as _bpy
+import os as _os
 
 
 def expandpath(path):
     if path.startswith("//"):
-        return os.path.join(os.path.dirname(bpy.data.filename), path[2:])
+        return _os.path.join(_os.path.dirname(_bpy.data.filename), path[2:])
 
     return path
 
@@ -47,67 +52,74 @@ _unclean_chars = ''.join([chr(i) for i in _unclean_chars])
 
 
 def clean_name(name, replace="_"):
-    '''
+    """
+    Returns a name with characters replaced that may cause problems under various circumstances, such as writing to a file.
     All characters besides A-Z/a-z, 0-9 are replaced with "_"
     or the replace argumet if defined.
-    '''
+    """
     for ch in _unclean_chars:
         name = name.replace(ch, replace)
     return name
 
 
 def display_name(name):
-    '''
-    Only capitalize all lowercase names, mixed case use them as is.
-    should work with filenames and module names.
-    '''
-    name_base = os.path.splitext(name)[0]
+    """
+    Creates a display string from name to be used menus and the user interface.
+    Capitalize the first letter in all lowercase names, mixed case names are kept as is.
+    Intended for use with filenames and module names.
+    """
+    name_base = _os.path.splitext(name)[0]
 
     # string replacements
     name_base = name_base.replace("_colon_", ":")
 
     name_base = name_base.replace("_", " ")
 
-    if name_base.lower() == name_base:
-        return ' '.join([w[0].upper() + w[1:] for w in name_base.split()])
+    if name_base.islower():
+        return name_base.capitalize()
     else:
         return name_base
 
 
 # base scripts
-_scripts = os.path.join(os.path.dirname(__file__), os.path.pardir, os.path.pardir)
-_scripts = (os.path.normpath(_scripts), )
+_scripts = _os.path.join(_os.path.dirname(__file__), _os.path.pardir, _os.path.pardir)
+_scripts = (_os.path.normpath(_scripts), )
 
 
 def script_paths(*args):
+    """
+    Returns a list of valid script paths from the home directory and user preferences.
+
+    Accepts any number of string arguments which are joined to make a path.
+    """
     scripts = list(_scripts)
 
     # add user scripts dir
-    user_script_path = bpy.context.user_preferences.filepaths.python_scripts_directory
+    user_script_path = _bpy.context.user_preferences.filepaths.python_scripts_directory
 
     if not user_script_path:
         # XXX - WIN32 needs checking, perhaps better call a blender internal function.
-        user_script_path = os.path.join(os.path.expanduser("~"), ".blender", "scripts")
+        user_script_path = _os.path.join(_os.path.expanduser("~"), ".blender", "scripts")
 
-    user_script_path = os.path.normpath(user_script_path)
+    user_script_path = _os.path.normpath(user_script_path)
 
-    if user_script_path not in scripts and os.path.isdir(user_script_path):
+    if user_script_path not in scripts and _os.path.isdir(user_script_path):
         scripts.append(user_script_path)
 
     if not args:
         return scripts
 
-    subdir = os.path.join(*args)
+    subdir = _os.path.join(*args)
     script_paths = []
     for path in scripts:
-        path_subdir = os.path.join(path, subdir)
-        if os.path.isdir(path_subdir):
+        path_subdir = _os.path.join(path, subdir)
+        if _os.path.isdir(path_subdir):
             script_paths.append(path_subdir)
 
     return script_paths
 
 
-_presets = os.path.join(_scripts[0], "presets") # FIXME - multiple paths
+_presets = _os.path.join(_scripts[0], "presets") # FIXME - multiple paths
 
 
 def preset_paths(subdir):
@@ -115,4 +127,4 @@ def preset_paths(subdir):
     Returns a list of paths for a spesific preset.
     '''
 
-    return (os.path.join(_presets, subdir), )
+    return (_os.path.join(_presets, subdir), )
index a1278f4f7cda89161ef846b81029e252c283cc6c..668951980e8ac3c3ed7b3fe5f301d9ad0512789c 100644 (file)
@@ -94,19 +94,19 @@ class _GenericBone:
     def x_axis(self):
         """ Vector pointing down the x-axis of the bone.
         """
-        return self.matrix.rotationPart() * Vector(1,0,0)
+        return self.matrix.rotation_part() * Vector(1,0,0)
     
     @property
     def y_axis(self):
         """ Vector pointing down the x-axis of the bone.
         """
-        return self.matrix.rotationPart() * Vector(0,1,0)
+        return self.matrix.rotation_part() * Vector(0,1,0)
     
     @property
     def z_axis(self):
         """ Vector pointing down the x-axis of the bone.
         """
-        return self.matrix.rotationPart() * Vector(0,0,1)
+        return self.matrix.rotation_part() * Vector(0,0,1)
 
     @property
     def basename(self):
@@ -236,7 +236,7 @@ class EditBone(StructRNA, _GenericBone):
         Expects a 4x4 or 3x3 matrix.
         """
         from Mathutils import Vector
-        z_vec = self.matrix.rotationPart() * Vector(0.0, 0.0, 1.0)
+        z_vec = self.matrix.rotation_part() * Vector(0.0, 0.0, 1.0)
         self.tail = matrix * self.tail
         self.head = matrix * self.head
         scalar = matrix.median_scale
@@ -421,8 +421,8 @@ class OrderedMeta(type):
 
 
 # Only defined so operators members can be used by accessing self.order
-class Operator(StructRNA, metaclass=OrderedMeta):
-    __slots__ = ()
+#class Operator(StructRNA, metaclass=OrderedMeta):
+#    __slots__ = ()
 
 
 class Macro(StructRNA, metaclass=OrderedMeta):
index 6fba34b7b58cd71477f8c08bda823aac20cc6e34..e10589d17d00adfac04f8c1c20f9c3f1e2b08631 100644 (file)
@@ -30,14 +30,14 @@ def get_hub(co, _hubs, EPS_SPLINE):
             if (hub.co - co).length < EPS_SPLINE:
                 return hub
 
-        key = co.toTuple(3)
+        key = co.to_tuple(3)
         hub = _hubs[key] = Hub(co, key, len(_hubs))
         return hub
     else:
         pass
 
         '''
-        key = co.toTuple(3)
+        key = co.to_tuple(3)
         try:
             return _hubs[key]
         except:
@@ -274,9 +274,7 @@ def get_splines(gp):
 
 
 def xsect_spline(sp_a, sp_b, _hubs):
-    from Mathutils import LineIntersect
-    from Mathutils import MidpointVecs
-    from Geometry import ClosestPointOnLine
+    from Geometry import ClosestPointOnLine, LineIntersect
     pt_a_prev = pt_b_prev = None
     EPS_SPLINE = (sp_a.length + sp_b.length) / (EPS_SPLINE_DIV * 2)
     pt_a_prev = sp_a.points[0]
@@ -296,7 +294,7 @@ def xsect_spline(sp_a, sp_b, _hubs):
                         # if f >= 0.0-EPS_SPLINE and f <= 1.0+EPS_SPLINE:
                         if f >= 0.0 and f <= 1.0:
                             # This wont happen often
-                            co = MidpointVecs(xsect[0], xsect[1])
+                            co = xsect[0].lerp(xsect[1], 0.5)
                             hub = get_hub(co, _hubs, EPS_SPLINE)
 
                             sp_a.hubs.append((a, hub))
@@